---
title: "Harvard STEPS CVD risk"
author: "Arpita Ghosh"
date: "18 Mar, 2020"
output: html_notebook
---

```{r load libraries}
rm(list=ls())
library(forestplot)
library(ggrepel)
library(gridExtra)
library(geepack)
library(lme4)
library(wesanderson)
library(foreign)
library(grid)
library(irr)
library(naniar)
library(tidyverse)
library(haven)
library(dplyr)
library(magrittr)
library(srvyr)
library(survey)
library(data.table)
library("openxlsx")
library(gtable)
library(forcats)
windowsFonts(Times=windowsFont("TT Times New Roman"))
```

```{r read in the data, select relevant variables, exclude observations}
# Import data #######################################################################
STEPS<-read_dta(file="DATASETS/HPACC-2019_10_16-CLEAN.dta")
dim(STEPS)
tabeadat<-read.csv(file="DATASETS/dat.analysis.tabea.csv",header=TRUE)
dim(tabeadat)
tabeadat <- tabeadat %>%   mutate(hh_id = ifelse(is.na(hh_id), p_id, hh_id))

#Select relevant variables 
data<-STEPS %>% dplyr::select(Country, year, svy, w2, sex, age, age_5yr, age_10yr, bmi, bmicat, csmoke, daily_smoke, past_smk,countryRegionclass, psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg, clin_hypt, bp_ms_new, hypt_new, hypt_med_new, 
ex_brazil_v0025,bp_ms,bpms_hypt_adv_med,wealth_quintile,edyears, educat,marital, working,fast_new,fbg_new,any_med,hba1c_p)
data$Country[data$Country=="Azerbaijan"]<-"Azerbaijan0"
data<-data[!(data$Country %in% c("Benin","Kyrgyz Republic")), ]

# There are 31 obs with a missing psu_num,  these get deleted
dim(data)
table(data$Country[is.na(data$psu_num)])
#data <- data %>%   filter(is.na(psu_num)==F)
data<-data[!is.na(data$psu_num),]
dim(data)

# Create a new individual participant ID that is truly unique
data <- data %>%   mutate(p_id = row_number())

# Restrict dataset in Brazil to those who were eligible for a BP measurement
data <- data %>%
  filter(ex_brazil_v0025==1 | is.na(ex_brazil_v0025==T))

# Restrict dataset in Ecuador to those who were eligible for the BP measurement (45k were skipped because of skip pattern)
data <- data %>% 
  mutate(ineligibleecuador = ifelse(Country=="Ecuador" & 
                                      (bp_ms == "999999999" | bp_ms == "666666666"), 1, 0)) %>% 
  filter(ineligibleecuador != 1)
data$educat <- as.factor(data$educat )
data <- data %>% 
  mutate(educat_labdetail = fct_recode(educat,"No formal schooling" = "0", "Some primary school" = "1","Primary school completed" = "2","Some high school" = "3","High school or above" ="4"))
data$edulte4yr<-with(data,ifelse(educat_labdetail %in% c("No formal schooling","Some primary school" ),1,0))

data$wpartner<-with(data,ifelse(marital ==2|marital ==6,1,ifelse(marital %in% c(1,3:5,7),0,NA)))

# final diabetes variable: 'diabetes' - this matches with clin_dia, except for India
#dibetes lancet paper - "Participants missing data on fasting status were assumed to be fasting as they were told to fast as part of all survey protocols, except in India, where participants were not instructed to fast"
data$fast_new_rev<-data$fast_new
data$fast_new_rev[data$Country=="India" & is.na(data$fast_new)]<-0
data$diabetes=0
data$diabetes[(data$fbg_new>=7.0 & data$fbg_new<800 & (data$fast_new_rev>=1 | is.na(data$fast_new_rev))) | (data$fbg_new>=11.1 & data$fbg_new<800 & data$fast_new_rev==0) | data$any_med==1]<-1
data$diabetes[((data$hba1c_p>=6.5 & data$hba1c_p<800) | data$any_med==1) & (data$Country %in% c("Fiji","Indonesia","Mexico","South Africa"))]<-1
data$diabetes[(data$fbg_new>800 | is.na(data$fbg_new>800)) & !(data$Country %in% c("Fiji","Indonesia","Mexico","South Africa"))]<-NA
data$diabetes[(data$hba1c_p>800 | is.na(data$hba1c_p>800))& (data$Country %in% c("Fiji","Indonesia","Mexico","South Africa"))]<-NA

data<-data %>% dplyr::select(Country, year, svy, w2, 
sex, age, age_5yr, age_10yr, 
        bmi, bmicat, csmoke, daily_smoke, past_smk,countryRegionclass, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_new,hypt_med_new,p_id,wealth_quintile,edyears,edulte4yr,wpartner,working,diabetes)
data<-data[data$Country!="South Africa DHS",]
oldcountries<-unique(data$Country)

```

```{r : add new countries}
cols<-c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")
newcountries<-c("Algeria","Azerbaijan","Belarus","Botswana","Cambodia","Moldova","Tajikstan","Kyrgyz Republic","Morocco","Sudan","Benin")
sapply(tabeadat[tabeadat$Country %in% newcountries,cols],class)
#table(droplevels(tabeadat[tabeadat$Country %in% newcountries,]$year),useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$year,useNA="ifany")
#table(droplevels(tabeadat[tabeadat$Country %in% newcountries,]$svy),useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$svy,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$sex,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$age,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$age_5yr,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$age_10yr,useNA="ifany")
#table(droplevels(tabeadat$Country[tabeadat$Country %in% newcountries & is.na(tabeadat$bmi)]))
table((tabeadat$Country[tabeadat$Country %in% newcountries & is.na(tabeadat$bmi)]))
range(tabeadat[tabeadat$Country %in% newcountries,]$bmi,na.rm=TRUE)
table(tabeadat[tabeadat$Country %in% newcountries,]$bmicat,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$csmoke,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$past_smk,useNA="ifany")
range(tabeadat[tabeadat$Country %in% newcountries,]$sbp_avg,na.rm=TRUE)
range(tabeadat[tabeadat$Country %in% newcountries,]$dbp_avg,na.rm=TRUE)
#table(droplevels(tabeadat$Country[tabeadat$Country %in% newcountries & is.na(tabeadat$sbp_avg)]))
table((tabeadat$Country[tabeadat$Country %in% newcountries & is.na(tabeadat$sbp_avg)]))
table((tabeadat$Country[tabeadat$Country %in% newcountries & is.na(tabeadat$dbp_avg)]))
table(tabeadat[tabeadat$Country %in% newcountries,]$clin_hypt,useNA="ifany")
table(tabeadat[tabeadat$Country %in% newcountries,]$bp_ms_new,useNA="ifany")
table((tabeadat$Country[tabeadat$Country %in% newcountries & tabeadat$bp_ms_new==2]))
table(tabeadat[tabeadat$Country %in% newcountries,]$hypt_med_new,useNA="ifany")

#Algeria
tempdat<-read_dta(file="DATASETS/Algeria 2016/dza2016.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco products, such as cigarettes, cigars or pipes?, Do you currently smoke tobacco products daily?  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
#Have you ever had your blood pressure measured by a doctor or other health worker? Have you ever been told by a doctor or other health worker that you have raised blood pressure or hypertension? 
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

table(tempdat$b1,is.na(tempdat$b5),useNA="ifany")
tempdat$fast<-with(tempdat,ifelse(b1==2 |is.na(b1),1,ifelse(b1==1,0,NA)))# b1 missing = fasting
tempdat$fast_NA<-with(tempdat,ifelse(is.na(b1),NA,ifelse(b1==2 ,1,ifelse(b1==1,0,NA))))# b1 missing = fasting missing
#>=7.0 mmol/l (126 mg/dl), random plasma glucose >=11.1 mmol/l (200 mg/dl)
range(tempdat$b5,na.rm=T)
tempdat$highbg<-with(tempdat,ifelse((fast==1 & b5>=126 ) | (fast==0 & b5>=200),1,0))
with(tempdat[tempdat$fast==1,],tapply(b5,highbg,range,na.rm=T));with(tempdat[tempdat$fast==0,],tapply(b5,highbg,range,na.rm=T))
tempdat$diab_med<-with(tempdat,ifelse(h6==1 & h7a==1 & (h8==1|h9==1),1,0))
tempdat$diabetes<-with(tempdat,ifelse(highbg==1 |diab_med==1,1,0))
table(tempdat$diabetes,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Algeria",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg","dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Algeria") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Algeria"
newdat$countryRegionclass<-1 # "Africa"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-newdat

#Azerbaijan
tempdat<-read_dta(file="DATASETS/Azerbaijan 2017/aze2017.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
range(tempdat$bmi,na.rm=T)
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Azerbaijan",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Azerbaijan") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Azerbaijan"
newdat$countryRegionclass<-4 # "EEM"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Belarus
tempdat<-read_dta(file="DATASETS/Belarus 2016/blr2016.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Belarus",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Belarus") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Belarus"
newdat$countryRegionclass<-4 # "EEM"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Benin 2015
tempdat<-read_dta(file="DATASETS/Benin 2015/ben2015.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Benin",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Benin") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Benin"
newdat$countryRegionclass<-1 # "Africa"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Botswana
tempdat<-read_dta(file="DATASETS/Botswana 2014/bwa2014.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:4,1,ifelse(c8 %in% 5:10,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

 
round(100*apply(is.na(tabeadat[tabeadat$Country=="Botswana",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Botswana") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Botswana"
newdat$countryRegionclass<-1 # "Africa"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Cambodia 2010
tempdat<-read_dta(file="DATASETS/Cambodia 2010/khm2010.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m4:weight(kg); m3: height(cm)
table(is.na(tempdat$m4));range(tempdat$m4,na.rm=T);head(table(tempdat$m4));tail(table(tempdat$m4))
table(is.na(tempdat$m3));range(tempdat$m3,na.rm=T);head(table(tempdat$m3));tail(table(tempdat$m3))
tempdat$wt_kg<-ifelse(tempdat$m4>600,NA,tempdat$m4)
tempdat$ht_cm<-ifelse(tempdat$m3>600,NA,tempdat$m3)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3a,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3a==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

 
round(100*apply(is.na(tabeadat[tabeadat$Country=="Cambodia",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Cambodia") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Cambodia"
newdat$countryRegionclass<-3 # "Southeast Asia and the western Pacific"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Kyrgyzstan 2013
tempdat<-read_dta(file="DATASETS/Kyrgyzstan 2013/kgz2013.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
range(tempdat$bmi,na.rm=T)
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Kyrgyz Republic",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Kyrgyz Republic") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Kyrgyzstan"
newdat$countryRegionclass<-4 # "EEM"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Morocco 2017
tempdat<-read_dta(file="DATASETS/Morocco 2017/mar2017.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
range(tempdat$bmi,na.rm=T)
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Morocco",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Morocco") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Morocco"
newdat$countryRegionclass<-1 # "Africa"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

#Moldova 2013
tempdat<-read_dta(file="DATASETS/Republic of Moldova 2013/Republic of Moldova.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
range(tempdat$bmi,na.rm=T)
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Moldova",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Moldova") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Moldova"
newdat$countryRegionclass<-4 # "EEM"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

##Sudan
tempdat<-read_dta(file="DATASETS/Sudan/sdn2016.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
range(tempdat$bmi,na.rm=T)
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))
tempdat$edyears_1<-ifelse(tempdat$c4==77,NA,tempdat$c4)
tempdat$edyears_1[tempdat$xc1==2]<-0

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Sudan",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Sudan") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Sudan"
newdat$countryRegionclass<-1 # "Africa"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
#newdat$edyears_1<-tempdat$edyears_1[match(newdat$p_id,tempdat$pid)]
#table(newdat$edyears_1,newdat$edyears,useNA="ifany")
newdat$edyears<-tempdat$edyears_1[match(newdat$p_id,tempdat$pid)]
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$edulte4yr[newdat$xc1==2]<-1
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

##Tajikistan 2016
tempdat<-read_dta(file="DATASETS/Tajikistan 2016/tjk2016.dta")
table(is.na(tempdat$pid));table(table(tempdat$pid))
table(tempdat$t1,tempdat$t2,useNA="ifany")#Do you currently smoke any tobacco  
tempdat$daily_smoke<-ifelse(tempdat$t1==2 | tempdat$t2==2,0,1)
table(tempdat$h1,tempdat$h2a,useNA="ifany")
tempdat$bp_ms_new<-with(tempdat,ifelse(h1==1,1,0))
tempdat$hypt_new<-with(tempdat,ifelse(h1==1 & h2a==1,1,0))
table(is.na(tempdat$wstep2))
#m12:weight(kg); m11: height(cm)
table(is.na(tempdat$m12));range(tempdat$m12,na.rm=T);head(table(tempdat$m12));tail(table(tempdat$m12))
table(is.na(tempdat$m11));range(tempdat$m11,na.rm=T);head(table(tempdat$m11));tail(table(tempdat$m11))
tempdat$wt_kg<-ifelse(tempdat$m12>600,NA,tempdat$m12)
tempdat$ht_cm<-ifelse(tempdat$m11>600,NA,tempdat$m11)
tempdat$bmi<-tempdat$wt_kg/(tempdat$ht_cm/100)^2
range(tempdat$bmi,na.rm=T)
tempdat$bmi[tempdat$bmi>80]<-NA
tempdat$bmicat<-with(tempdat,ifelse(bmi < 18.5, 1,ifelse(bmi >= 18.5 & bmi < 25, 2,ifelse(bmi >= 25 & bmi <30, 3,ifelse(bmi >= 30 & bmi<=80, 4, NA)))))
with(tempdat,table(h1,h2a,h3,useNA="ifany"))
tempdat$hypt_med_new_1<-with(tempdat,ifelse(h1==1 & h2a==1 & h3==1,1,0))

table(tempdat$c8,useNA="ifany")
print_labels(tempdat$c8, name = NULL)
tempdat$working<-with(tempdat,ifelse(c8 %in% 1:3,1,ifelse(c8 %in% 4:9,0,NA)))
table(tempdat$c8,tempdat$working,useNA="ifany")

table(tempdat$c7,useNA="ifany")
print_labels(tempdat$c7, name = NULL)
tempdat$wpartner<-with(tempdat,ifelse(c7 ==2|c7 ==6,1,ifelse(c7 %in% c(1,3:5),0,NA)))
table(tempdat$c7,tempdat$wpartner,useNA="ifany")

round(100*apply(is.na(tabeadat[tabeadat$Country=="Tajikistan",c("year", "svy", "sex", "age", "age_5yr", "age_10yr","bmi", "bmicat", "csmoke", "past_smk","psu_num", "hh_id", "stratum_num","Population2015","sbp_avg", "dbp_avg","clin_hypt","bp_ms_new","hypt_med_new","p_id","wealth_quintile","edyears")]),2,mean),1)
newdat<-tabeadat %>% filter(Country=="Tajikistan") %>% dplyr::select(year, svy, 
sex, age, age_5yr, age_10yr,bmi, bmicat, csmoke, past_smk, 
        psu_num, hh_id, stratum_num, Population2015,sbp_avg, dbp_avg,
        clin_hypt,bp_ms_new,hypt_med_new,p_id,wealth_quintile,edyears)
table(newdat$bp_ms_new,tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)],useNA="ifany");table(newdat$bmicat,tempdat$bmicat[match(newdat$p_id,tempdat$pid)],useNA="ifany");range(newdat$bmi,na.rm=T);range(tempdat$bmi[match(newdat$p_id,tempdat$pid)],na.rm=T)
newdat$Country<-"Tajikistan"
newdat$countryRegionclass<-4 # "EEM"
table(is.na(match(newdat$p_id,tempdat$pid)))
newdat$w2<-tempdat$wstep2[match(newdat$p_id,tempdat$pid)]
newdat$daily_smoke<- tempdat$daily_smoke[match(newdat$p_id,tempdat$pid)]
newdat$bp_ms_new<-tempdat$bp_ms_new[match(newdat$p_id,tempdat$pid)]
newdat$hypt_new<-tempdat$hypt_new[match(newdat$p_id,tempdat$pid)]
newdat$bmi<-tempdat$bmi[match(newdat$p_id,tempdat$pid)]
newdat$bmicat<-tempdat$bmicat[match(newdat$p_id,tempdat$pid)]
newdat$hypt_med_new_1<-tempdat$hypt_med_new_1[match(newdat$p_id,tempdat$pid)]
table(newdat$hypt_med_new,newdat$hypt_med_new_1,useNA="ifany")
newdat$working<-tempdat$working[match(newdat$p_id,tempdat$pid)]
newdat$wpartner<-tempdat$wpartner[match(newdat$p_id,tempdat$pid)]
#newdat$diabetes<-tempdat$diabetes[match(newdat$p_id,tempdat$pid)]
newdat$edulte4yr<-with(newdat,ifelse(edyears<5,1,0))
newdat$p_id<-paste(newdat$Country,newdat$p_id,sep="_")
newdata<-rbind(newdata,newdat)

round(100*apply(is.na(newdata),2,mean),1)
table(newdata$Country,newdata$countryRegionclass,useNA="ifany")
#table(droplevels(newdata$year),useNA="ifany")
table(newdata$year,useNA="ifany")
#table(droplevels(newdata$svy),useNA="ifany")
table(newdata$svy,useNA="ifany")
table(newdata$sex,useNA="ifany")
table(newdata$age,useNA="ifany")
table(newdata$age_5yr,useNA="ifany")
table(newdata$age_10yr,useNA="ifany")
table(newdata$Country[is.na(newdata$bmi)])
range(newdata$bmi,na.rm=TRUE)
table(newdata$bmicat,useNA="ifany")
table(newdata$csmoke,useNA="ifany")
table(newdata$past_smk,useNA="ifany")
table(newdata$daily_smoke,useNA="ifany")
range(newdata$sbp_avg,na.rm=TRUE)
range(newdata$dbp_avg,na.rm=TRUE)
table(newdata$Country[is.na(newdata$sbp_avg)])
table(newdata$Country[is.na(newdata$dbp_avg)])
table(newdata$clin_hypt,useNA="ifany")
table(newdata$bp_ms_new,useNA="ifany")
table(newdata$hypt_new,useNA="ifany")
table(newdata$hypt_med_new,useNA="ifany")
table(newdata$hypt_med_new_1,useNA="ifany")
table(newdata$edulte4yr,useNA="ifany")
table(newdata$wealth_quintile,useNA="ifany")

newdata$hypt_med_new<-newdata$hypt_med_new_1
dim(newdata)
newdata<-newdata %>% dplyr::select(-hypt_med_new_1)
dim(newdata)
newdata$diabetes<-NA
data<-rbind(data,newdata)

```

```{r : manage data}
#data<-STEPS
# Rename Timor-Leste
#data <- data %>%   mutate(Country = ifelse(Country=="Timor Leste", "Timor-Leste", Country))
data$Country[data$Country=="Timor Leste"]<-"Timor-Leste"
data$Country[data$Country=="St. Vincent & the Grenadines"]<-"StVG"

data$ccodes<- with(data,ifelse(Country=="Ukraine","UKR",ifelse(Country=="South Africa","ZAF",ifelse(Country=="Seychelles","SYC",ifelse(Country=="Romania","ROU",ifelse(Country=="Peru","PER",ifelse(Country=="Fiji","FJI",ifelse(Country=="Egypt","EGY",ifelse(Country=="Bangladesh","BGD",ifelse(Country=="Azerbaijan","AZE",(ifelse(Country=="Albania","ALB",ifelse(Country=="Belize","BLZ",ifelse(Country=="Benin","BEN",ifelse(Country=="Bhutan","BTN",ifelse(Country=="Brazil","BRA",ifelse(Country=="Burkina Faso","BFA",ifelse(Country=="Chile","CHL",ifelse(Country=="China","CHN",ifelse(Country=="Comoros","COM",ifelse(Country=="Costa Rica","CRI",ifelse(Country=="Ecuador","ECU",ifelse(Country=="Georgia","GEO",ifelse(Country=="Ghana","GHA",ifelse(Country=="Grenada","GRD",ifelse(Country=="Guyana","GUY",ifelse(Country=="India","IND",ifelse(Country=="Indonesia","IDN",ifelse(Country=="Kazakhstan","KAZ",ifelse(Country=="Kenya","KEN",ifelse(Country=="Kyrgyz Republic","KGZ",ifelse(Country=="Lebanon","LBN",ifelse(Country=="Lesotho","LSO",ifelse(Country=="Liberia","LBR",ifelse(Country=="Mexico","MEX",ifelse(Country=="Mongolia","MNG",ifelse(Country=="Mozambique","MOZ",ifelse(Country=="Namibia","NAM",ifelse(Country=="Nepal","NPL",ifelse(Country=="Russian Federation","RUS",ifelse(Country=="StVG","VCT",ifelse(Country=="Swaziland","SWZ",ifelse(Country=="Tanzania","TZA",ifelse(Country=="Timor-Leste","TLS",ifelse(Country=="Togo","TGO",ifelse(Country=="Uganda","UGA",ifelse(Country=="Zanzibar","TZA",ifelse(Country=="Algeria","DZA",ifelse(Country=="Azerbaijan","AZE",ifelse(Country=="Belarus","BLR",NA))))))))))))))))))))))))))))))))))))))))))))))))))
data$ccodes[data$Country=="Benin2015"]<-"BEN"
data$ccodes[data$Country=="Botswana"]<-"BWA";data$ccodes[data$Country=="Cambodia"]<-"KHM";data$ccodes[data$Country=="Kyrgyzstan"]<-"KGZ";data$ccodes[data$Country=="Morocco"]<-"MAR";data$ccodes[data$Country=="Sudan"]<-"SDN";data$ccodes[data$Country=="Tajikistan"]<-"TJK";data$ccodes[data$Country=="Moldova"]<-"MDA"
data$ccodes[data$Country=="Vanuatu"]<-"VUT"

# revise sbp var
data <- dplyr::mutate(data, sbp_avg_rev = ifelse(sbp_avg>=999, NA, sbp_avg)) # assign these values as NA: 666666688, 777777792, 888888896, 1e+09
# revise dbp var
data <- dplyr::mutate(data, dbp_avg_rev = ifelse(dbp_avg>=999, NA, dbp_avg)) # assign these values as NA: 666666688, 777777792, 888888896, 1e+09


data$controlled_hypt<-NA
data$controlled_hypt[data$sbp_avg_rev<140 & data$dbp_avg_rev<90]<-1
data$controlled_hypt[(data$sbp_avg_rev>=140 & data$dbp_avg_rev>=90) | (data$sbp_avg_rev>=140 & is.na(data$dbp_avg_rev)) | (is.na(data$sbp_avg_rev) & data$dbp_avg_rev>=90) | (is.na(data$sbp_avg_rev) & is.na(data$dbp_avg_rev))]<-0

data$bpms_hypt_med<-0
data$bpms_hypt_med[data$bp_ms_new==1 & data$hypt_new==1 & data$hypt_med_new==1 ]<-1
data$bpms_hypt_med[data$hypt_new==1 & data$hypt_med_new==1 & data$Country %in% c("Albania", "China", "Egypt", "Ghana", "Indonesia","Kazakhstan", "Kyrgyz Republic","Mexico","Peru","Russian Federation","Seychelles")]<-1

data$bpms_hypt_med_ctl<-0
data$bpms_hypt_med_ctl[data$bp_ms_new==1 & data$hypt_new==1 & data$hypt_med_new==1 & data$controlled_hypt==1]<-1
data$bpms_hypt_med_ctl[data$hypt_new==1 & data$hypt_med_new==1 & data$controlled_hypt==1 & (data$Country=="Albania" | data$Country=="China" | data$Country=="Egypt" | data$Country=="Ghana" | data$Country=="Indonesia" | data$Country=="Kazakhstan" 	| data$Country=="Kyrgyz Republic" | data$Country=="Mexico" | data$Country=="Peru" | data$Country=="Russian Federation" | data$Country=="Seychelles")]<-1

# set relevant vars to missing for countries in which some steps are missing
data <- data %>% 
  mutate(bpms_hypt_med_ctl = ifelse(Country=="Romania" | Country=="Seychelles", NA, bpms_hypt_med_ctl),bpms_hypt_med = ifelse(Country=="Romania" | Country=="Seychelles", NA, bpms_hypt_med))

#Create new variables
#data$Country <- as.factor(data$Country)

data <- data %>% 
  mutate(region=
    ifelse(countryRegionclass ==1, "Africa", 
        ifelse(countryRegionclass ==2, "Latin America and the Caribbean",
            ifelse(countryRegionclass ==3, "Southeast Asia and the western Pacific",
                ifelse(countryRegionclass ==4, "Europe and the eastern Mediterranean", NA)))))
data$region[data$region=="Europe and the eastern Mediterranean"]<-"EEM"
data$region<-factor(data$region,levels=c("Latin America and the Caribbean",   
"EEM", "Southeast Asia and the western Pacific",
"Africa"))

# Make sure psu_num really is unique within each country
data <- data %>%   mutate(psu_num = str_c(Country, psu_num, sep = "_")) 

# clean hh_id variable
data <- data %>%   mutate(hh_id = ifelse(hh_id=="666666666" | hh_id == "666666688", p_id, hh_id)) # these are the only two missing/weird codes in that var; missing hh_id is here set to individual participant id


# Add labels to main predictors
data$bmicat <- as.factor(data$bmicat )

data <- data %>% 
  mutate(female_lab = as.factor(
                          ifelse(sex=="1", "Female", 
                             ifelse(sex=="0", "Male", NA))),
         
         age_10yr_lab = as.factor(
                     ifelse(age_10yr=="1", "15-24",
                        ifelse(age_10yr=="2", "25-34",
                            ifelse(age_10yr=="3", "35-44",
                                ifelse(age_10yr=="4", "45-54",
                                     ifelse(age_10yr=="5", "55-64",
                                          ifelse(age_10yr=="6", "65+", NA
                                                 ))))))),

         bmicat_lab = fct_recode(bmicat,
                        "<18.5" = "1", 
                        "18.5-24.9" = "2",
                        "25.0-29.9" = "3",
                        ">=30.0" = "4")
)

# Framingham risk score equation applies to only 30-74 year olds
data$age30to74=0
data$age30to74[data$age >=30 & data$age<75]<-1 # age missing assigned 0
#data$age30to74[data$age >=30 & data$age<75 & (data$pregnant!=1 | is.na(data$pregnant))]<-1

#WHO2019 risk charts apply to 40 to 74 years
data$age40to74=0
data$age40to74[data$age >=40 & data$age<75]<-1 # age missing assigned 0

#The estimated effect of one drug at standard dose in lowering blood pressure from a pretreatment blood pressure P is therefore [9.1+0.10(P-154)] systolic and [5.5+0.11(P???97)] diastolic
#postbp=prebp-9.1-0.10(prebp-154)
#p'=9.1+0.10(P-154)

data$sbp_adjtrt<-data$sbp_avg_rev
data$sbp_adjtrt[data$bpms_hypt_med==1 & !is.na(data$bpms_hypt_med) & !is.na(data$sbp_avg_rev)]<-(data$sbp_avg_rev[data$bpms_hypt_med==1 & !is.na(data$bpms_hypt_med) &  !is.na(data$sbp_avg_rev)]+9.1-0.1*154)/(1-0.1)
data$sbp_adjtrt[is.na(data$bpms_hypt_med)]<-NA

data$dbp_adjtrt<-data$dbp_avg_rev
data$dbp_adjtrt[data$bpms_hypt_med==1 & !is.na(data$bpms_hypt_med) & !is.na(data$dbp_avg_rev)]<-(data$dbp_avg_rev[data$bpms_hypt_med==1 & !is.na(data$bpms_hypt_med) & !is.na(data$dbp_avg_rev)]+5.5-0.11*97)/(1-0.11)
data$dbp_adjtrt[is.na(data$bpms_hypt_med)]<-NA


#'smoke' from Pascal's code: htncascades_2018-11-11.Rmd
data$csmoke[data$Country=="Indonesia" & data$csmoke=="999999999"]<-0
data$past_smk[data$Country=="Chile" & data$past_smk=="999999999"]<-1

data <- data %>% 
  mutate(past_smk_nomiss = ifelse(is.na(past_smk)==T, 989898, past_smk),
         csmoke_nomiss = ifelse(is.na(csmoke)==T, 989898, csmoke),
         psmk_P = ifelse(csmoke_nomiss=="1" & past_smk_nomiss !="555555584" & past_smk_nomiss !="666666666" & past_smk_nomiss !="666666688", 1, past_smk_nomiss),
         csmk_P = ifelse(psmk_P==0 & Country=="Ecuador", 0, csmoke_nomiss), 
         smoke_detail = ifelse(csmk_P==1, "current smoker",
                ifelse(psmk_P==1 & csmk_P==0, "past smoker",
                      ifelse(csmk_P<2 | psmk_P<2, "never smoked", NA))),
         smoke =ifelse(is.na(smoke_detail), NA, 
                ifelse(smoke_detail=="current smoker", 1, 0)))

#'csmk' from Michaela's code
smokdat<-data[,c("past_smk","csmoke","daily_smoke","Country")]

with(smokdat[smokdat$Country=="China",],table(csmoke,past_smk,useNA="ifany"))
smokdat$csmoke[smokdat$Country=="China" & smokdat$past_smk==0 & smokdat$csmoke == 999999999 & !is.na(smokdat$csmoke) & !is.na(smokdat$past_smk)]<-0  #China : 7106 values for which past_smk=0 and csmoke = 999999999 (i.e. skipped), is coded as 0

smokdat$past_smk[smokdat$past_smk %in% c(666666666,666666688,555555584)]<-".m"
smokdat$past_smk[smokdat$Country=="India"]<-".m"
smokdat$past_smk[smokdat$Country %in% names(table(smokdat$Country))[tapply((smokdat$past_smk==1 |smokdat$past_smk==0)  & !is.na(smokdat$past_smk),smokdat$Country,sum)==0]] <-".m"
smokdat$csmk = smokdat$csmoke
smokdat$csmk[smokdat$Country=="Indonesia"]<- with(smokdat[smokdat$Country=="Indonesia",],ifelse(csmoke==999999999,0,csmoke)) # was coded wrongly in cleaning process 
smokdat$psmk<-with(smokdat,ifelse(csmk == 1 & !is.na(csmk == 1) & (past_smk!=".m" | is.na(past_smk)) ,1,past_smk))
smokdat$psmk[smokdat$Country=="Indonesia"]=with(smokdat[smokdat$Country=="Indonesia",],ifelse(past_smk == 999999999 | past_smk == 1e+09,0,psmk)) #was coded wrongly in cleaning process 
  
smokdat$psmk[smokdat$Country=="Chile"]<-with(smokdat[smokdat$Country=="Chile",],ifelse(past_smk == 999999999 | past_smk == 1e+09,1,psmk)) #in the cleaning code of the Chile raw data, observations who currently smoke where coded as ".s". We correct this here
smokdat$csmk[smokdat$Country=="Ecuador"]<-with(smokdat[smokdat$Country=="Ecuador",],ifelse(psmk == 0,0,csmk))#they have a different skip pattern

smokdat$daily_smoke[smokdat$daily_smoke %in% c(666666666,666666688,555555584)]<-".m"
smokdat$daily_smoke[smokdat$Country=="India"]<-".m"
smokdat$daily_smoke[smokdat$Country %in% names(table(smokdat$Country))[tapply((smokdat$daily_smoke==1 |smokdat$daily_smoke==0)  & !is.na(smokdat$daily_smoke),smokdat$Country,sum)==0]] <-".m"

smokdat$cdsmk = smokdat$daily_smoke
smokdat$cdsmk[smokdat$csmk == 0 & (smokdat$cdsmk != ".m"|is.na(smokdat$cdsmk))] <-0
smokdat$cdsmk [is.na(smokdat$csmk) & (smokdat$cdsmk == 999999999 | smokdat$cdsmk == 1e+09)]<-NA 
smokdat$cdsmk[(smokdat$cdsmk == 999999999 | smokdat$cdsmk == 1e+09) & smokdat$csmoke == 1] <-NA
smokdat$csmk [smokdat$cdsmk == 1]<- 1 
data$csmk<-smokdat$csmk
data$currsmok<-with(data,ifelse((csmk==0 | csmk==1 ) & !is.na(csmk),csmk,NA))
table(data$smoke,data$currsmok,useNA="ifany")

#complete cases
data<-mutate(data,complete=ifelse(!is.na(age) & !is.na(sex) & !is.na(currsmok) & !is.na(sbp_adjtrt) & !is.na(bpms_hypt_med) & !is.na(bmi),1,0))

#weights
data$CVDrisk<-NA
data$CVDrisk=with(data,ifelse(complete==1,1,NA))

#generate CVD risk weight variable
data$w_CVDrisk<-NA
data$w_CVDrisk[data$age30to74==1]=data$w2[data$age30to74==1] # restrict to adults 30-74 years
data$w_CVDrisk[is.na(data$CVDrisk)]=NA 
		
#generate average weight in country and replace missing weights with this
data<-data %>% group_by(Country) %>% mutate(mean_w_CVDrisk=mean(w_CVDrisk,na.rm=T)) 
data$w_CVDrisk[data$CVDrisk==1 & !is.na(data$CVDrisk) & is.na(data$w_CVDrisk) & data$age30to74==1 ]<-data$mean_w_CVDrisk[data$CVDrisk==1 & !is.na(data$CVDrisk) & is.na(data$w_CVDrisk) & data$age30to74==1 ]

#every country is weighted equally
data<-data %>% group_by(Country) %>% mutate(CVDrisk_obs=ifelse(!is.na(w_CVDrisk),sum(w_CVDrisk,na.rm=T),NA))
data$weq_CVDrisk=NA
data$weq_CVDrisk=with(data,ifelse(!is.na(w_CVDrisk),w_CVDrisk/CVDrisk_obs,NA))

#sum weights by country and rescale weights by population size
sumpop2015<-sum(data$Population2015[!duplicated(data$Country) & !(data$Country %in% names(table(data$Country))[tapply(data$weq_CVDrisk,data$Country,sum,na.rm=T)==0] )])
data$wpop2015_CVDrisk=data$weq_CVDrisk*data$Population2015/sumpop2015

# world population
# data available at http://ghdx.healthdata.org/record/ihme-data/gbd-2017-population-estimates-1950-2017

popdat<-read.csv("DATASETS/IHME_GBD_2017_POP_2015_2017_Y2018M11D08.csv",header = TRUE)

#alternate weights
data$weights<-data$w2
data<-data %>% group_by(Country) %>% mutate(nobs=ifelse(!is.na(weights),sum(weights,na.rm=T),NA))
data$weights_eq=NA
data$weights_eq=with(data,ifelse(!is.na(weights),weights/nobs,NA))

data <- data %>% 
  mutate(age_WHO2019=floor(age),
         age_WHO2019_trunc=ifelse(age<40,40,age_WHO2019),
         gdr_WHO2019=ifelse(sex==1,2,ifelse(sex==0,1,NA)),
         smk_WHO2019=ifelse(currsmok==1,1,ifelse(currsmok==0,0,NA)),
         sbp_WHO2019=sbp_avg_rev,
         sbpadj_WHO2019=sbp_adjtrt,
         bmi_WHO2019=bmi
         )
```

```{r : create dat & missdat- Figure S3 }
#drop particiapnts aged < 30 years and > 74 years
dat<-filter (data,age30to74==1)
dim(dat)
#fdat<-filter(dat,sex==1)
#mdat<-filter(dat,sex==0)

##########   missing data
missdat<-dat %>% group_by(Country) %>% summarise_at(.vars= vars( sex, currsmok,bmi,sbp_adjtrt), funs(miss=sum(is.na(.))))%>% left_join(.,dat%>% group_by(Country) %>% summarise(n=n(),n.complete=sum(!is.na(sex) & !is.na(currsmok) & !is.na(bmi) & !is.na(sbp_adjtrt)  ))) %>% mutate_at(.vars= vars( sex_miss, currsmok_miss,bmi_miss,sbp_adjtrt_miss,n.complete),funs(paste(.," (",sprintf("%.1f",.*100/n),")",sep=""))) 

p<-with(dat,gg_miss_fct(x = dat[,c("Country","sex", "currsmok","bmi","sbp_adjtrt")], fct = Country))
p<-p+scale_y_discrete(limits=c("sex", "currsmok","bmi","sbp_adjtrt"),labels=str_wrap(c("Sex", "Smoking status","BMI","Systolic blood pressure","Taking antihypertensive drugs"),width=15))+labs(x="",y="",color="% Missing")+theme(text = element_text(size=15),axis.text = element_text(colour = "black",face = "bold"),plot.margin = margin(0,0.05 , -0.5, -0.5, "cm"))+coord_flip()+guides(fill=guide_legend(title="% Missing"))

#drop countries with 100% missing for one variable
dropCoun<-dat%>% group_by(Country) %>% summarise(n.complete=sum(!is.na(sex) & !is.na(currsmok) & !is.na(bmi) & !is.na(sbp_adjtrt) )) %>% filter(n.complete==0) %>% pull(Country)

dat1<-dat %>% filter(!(Country %in% dropCoun )) 
dat1<-dat1%>% filter(!is.na(w2) | complete==1)

dat1$Country<-with(dat1,ifelse(Country=="St. Vincent & the Grenadines","StVG",Country))
f <- function(x) factor(x, levels = unique(x))
dat1<-with(dat1, dat1[order(region, Country),]) 
dat1$Country<-f(dat1$Country)
p<-gg_miss_fct(x = dat1[,c("Country","sex", "currsmok","bmi","sbp_adjtrt")], fct = Country)
p<-p+scale_y_discrete(limits=c("sex", "currsmok","bmi","sbp_adjtrt"),labels=str_wrap(c("Sex", "Smoking status","BMI","Systolic blood pressure"),width=15))+labs(x="",y="",color="% Missing")+theme(text = element_text(size=16),axis.text = element_text(colour = "black",face = "bold"),plot.margin = margin(0,0.05 , -0.5, -0.5, "cm"),axis.text.x = element_text(
      angle = 90,
      hjust = 1,
      vjust = 0.5
  ))+guides(fill=guide_colorbar(title="% Missing"))+scale_fill_distiller(palette = "RdYlGn")

#scale_fill_distiller(palette = "Spectral")
#scale_fill_distiller(palette = "RdYlGn")

p$data$region<-dat1$region[match(p$data$Country,dat1$Country)]
p1<-p+facet_grid(.~ region, space = 'free_x', scales = 'free_x', switch = 'x',labeller = label_wrap_gen(width=17)) +  theme(strip.placement = 'outside',strip.background.x = element_blank(),strip.text = element_text(colour = "black",face = "bold",size=15)) 
p1
##################################################   complete cases
dat<-filter (dat,complete==1 )
#dat<-filter (dat,complete==1 & Country!="Zanzibar")
dim(dat)
```

```{r : Table 1 - sample characteristics }

table1.svy <- dat %>% 
            as_survey_design(stratum = stratum_num,
                        ids = psu_num, 
                        weights = weq_CVDrisk,
                        variables = c(region,Country,sex,female_lab,age,currsmok,bpms_hypt_med,sbp_avg_rev,bmi))

table1.wtd<-left_join(dat %>% group_by(Country,female_lab) %>% summarise(region=unique(region),n=n(),meanage=mean(age),medage=floor(median(age)),age_min=floor(min(age)),age_max=floor(max(age))) %>% mutate(agerange=paste(age_min,"-",age_max,sep="")) %>% select( region,Country,female_lab,n, agerange) %>% arrange(.,region,Country,female_lab),table1.svy %>% group_by(region,Country,female_lab) %>% summarise(meanage=survey_mean(age),medage=survey_quantile(age,0.5),pctsmok=survey_mean(currsmok, proportion=TRUE, prop_method="asin"),pctHTNmed=survey_mean(bpms_hypt_med, proportion=TRUE, prop_method="asin")) %>% mutate(meanage=sprintf("%.1f",meanage),medage_q50=floor(medage_q50),pctsmok=sprintf("%.1f",100*pctsmok),pctHTNmed=sprintf("%.1f",100*pctHTNmed)) %>% select(region,Country,female_lab,meanage,medage_q50,pctsmok,pctHTNmed) %>% arrange(.,region,Country,female_lab))
table1.wtd<-left_join(table1.wtd,table1.svy %>% subset(.,bpms_hypt_med==0) %>% group_by(region,Country,female_lab) %>% summarise(sbp_notrt=survey_quantile(sbp_avg_rev,c(0.25,0.5,0.75))) %>%  mutate (sbp_notrt=paste(sprintf("%.1f",sbp_notrt_q50)," (",sprintf("%.1f",sbp_notrt_q25),", ",sprintf("%.1f",sbp_notrt_q75),")",sep="")) %>% select(region,Country,female_lab,sbp_notrt) %>% arrange(.,region,Country,female_lab))
table1.wtd<-left_join(table1.wtd,table1.svy %>% subset(.,bpms_hypt_med==1) %>% group_by(region,Country,female_lab) %>% summarise(sbp_trt=survey_quantile(sbp_avg_rev,c(0.25,0.5,0.75))) %>%  mutate (sbp_trt=paste(sprintf("%.1f",sbp_trt_q50)," (",sprintf("%.1f",sbp_trt_q25),", ",sprintf("%.1f",sbp_trt_q75),")",sep="")) %>% select(region,Country,female_lab,sbp_trt) %>% arrange(.,region,Country,female_lab))
table1.wtd<-left_join(table1.wtd,table1.svy %>% group_by(region,Country,female_lab) %>% summarise(bmi=survey_quantile(bmi,c(0.25,0.5,0.75))) %>%  mutate (bmi=paste(sprintf("%.1f",bmi_q50)," (",sprintf("%.1f",bmi_q25),", ",sprintf("%.1f",bmi_q75),")",sep="")) %>% select(region,Country,female_lab,bmi) %>% arrange(.,region,Country,female_lab)) 
table1.wtd$region<-ifelse(!duplicated(table1.wtd$region),as.character(table1.wtd$region),"")
table1.wtd$Country<-with(table1.wtd,ifelse(!duplicated(Country),as.character(Country),""))
```

```{r : WHO2019}
df1<-read_dta(file="DATASETS/mydata_HPACC_28Apr2020_wrisk_40to74.dta")
dim(df1)
#df1$p_id<-as.character(df1$p_id)
#df1$p_id<-format(df1$p_id, scientific = FALSE, trim = TRUE)
df2<-read_dta(file="DATASETS/mydata_new_28Apr2020_wrisk_40to74.dta")
dim(df2)
WHO2019dat_wrisk<-rbind(df1,df2) 
dim(WHO2019dat_wrisk)
dat$WHO2019risk_40to74<-WHO2019dat_wrisk$cal2_who_cvdx_m2[match(dat$p_id,WHO2019dat_wrisk$p_id)]
dat$gbdreg<-WHO2019dat_wrisk$gbdreg[match(dat$p_id,WHO2019dat_wrisk$p_id)]
table(table(dat$p_id))
dat$pid<-1:dim(dat)[1]
```


```{r survey design}
dat$stratum_num <- as.factor(dat$stratum_num) 
dat$psu_num <- as.factor(dat$psu_num) 
dat$hh_id <- as.factor(dat$hh_id) 
dat$weq_CVDrisk <- as.numeric(dat$weq_CVDrisk) 
dat$wpop2015_CVDrisk <- as.numeric(dat$wpop2015_CVDrisk) 
dat$Country <- as.factor(dat$Country) 
dat$WHOrisk <- 100*as.numeric(dat$WHO2019risk_40to74) # can also be dat$WHO2019risk_30to74: same equation applied to 30-39 year olds
dat <- dat %>% 
  mutate(WHOriskcat = ifelse(WHOrisk<20, "Low", ifelse(WHOrisk<30,"Medium","High")))
dat <- dat %>% 
  mutate(WHOrisk5cat = ifelse(WHOrisk<5, "<5", ifelse(WHOrisk<10,"5-<10",ifelse(WHOrisk<20,"10-<20",ifelse(WHOrisk<30,"20-<30",">=30")))))
dat <- dat %>% 
  mutate(WHOrisk_low = ifelse(WHOrisk<20, 1, 0),WHOrisk_med = ifelse(WHOrisk>=20 & WHOrisk<30, 1, 0),WHOrisk_high = ifelse(WHOrisk>=30, 1, 0))
# Set psu_num in those countries with only one PSU equal to hh_id
onlyonepsu <- c("China", "Grenada", "Kazakhstan", "Mexico", "Romania", "Seychelles", "Uganda")
dat <- dat %>% 
  mutate(psu_num = ifelse(Country %in% onlyonepsu, hh_id, psu_num))

######## response to reviewers' comment:

table1.svy <- dat %>% 
            as_survey_design(stratum = stratum_num,
                        ids = psu_num, 
                        weights = weq_CVDrisk,
                        variables = c(region,Country,sex,female_lab,age,currsmok,bpms_hypt_med,sbp_avg_rev,bmi,diabetes))


tab1<-table1.svy %>% filter(table1.svy$variables$Country %in% oldcountries) %>% group_by(region,Country,female_lab) %>% summarise(meanage=survey_mean(age),pctsmok=survey_mean(currsmok, proportion=TRUE, prop_method="asin"),pctdiab=survey_mean(diabetes, proportion=TRUE, prop_method="asin",na.rm=T),pctHTNmed=survey_mean(bpms_hypt_med, proportion=TRUE, prop_method="asin")) %>% mutate(meanage=sprintf("%.1f",meanage),pctsmok=sprintf("%.1f",100*pctsmok),pctdiab=sprintf("%.1f",100*pctdiab),pctHTNmed=sprintf("%.1f",100*pctHTNmed)) %>% select(region,Country,female_lab,meanage,pctsmok,pctdiab,pctHTNmed) %>% arrange(.,region,Country,female_lab)
```

```{r : Figure 1  - WHOrisk distribution}
svy.WHOrisk <- dat %>% 
            filter(is.na(WHOrisk)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk,
                        variables = c(WHOrisk,WHOrisk5cat,WHOriskcat,
                                      Country, countryRegionclass, region,age, currsmok,bpms_hypt_med,sbp_avg_rev,bmi,bpms_hypt_med,
                                      sex, age_10yr,
                                      bmicat))
dat_WHOrisk <-svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% 
  summarize(meanWHOrisk_1 = survey_mean(WHOriskcat=="Low", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(WHOriskcat=="Medium", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_3 = survey_mean(WHOriskcat=="High", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se","meanWHOrisk_3_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_3),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_3)

dat_WHOrisk <- svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% 
  summarize(meanWHOrisk_1 = survey_mean(WHOrisk5cat=="<5", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(WHOrisk5cat=="5-<10", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_3 = survey_mean(WHOrisk5cat=="10-<20", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_4 = survey_mean(WHOrisk5cat=="20-<30", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_5 = survey_mean(WHOrisk5cat==">=30", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se","meanWHOrisk_3_se","meanWHOrisk_4_se","meanWHOrisk_5_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_5),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_5)



f <- function(x) factor(x, levels = unique(x))
orderCountry<-dat_WHOrisk[dat_WHOrisk$RiskCategory=="meanWHOrisk_1" & dat_WHOrisk$sex==0,] %>% arrange(desc(region),desc(Pct)) %>% select(Country) %>%  .$Country %>% f(.)
dat_WHOrisk<-dat_WHOrisk %>%
   arrange(desc(region), factor(Country, levels = levels(orderCountry)))


p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_1", "meanWHOrisk_2", "meanWHOrisk_3", "meanWHOrisk_4", "meanWHOrisk_5"),labels=c("<5%","5 - <10%","10 - <20%","20 - <30%",">=30%"))))) + geom_bar(stat="identity",width = 0.9)+labs(y="Percent",x="")+coord_flip()+guides(fill=guide_legend(title="10-year \nCVD Risk"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=10))+ scale_fill_manual(values=c("#990000","orange","yellow","darkgreen","green"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=15),axis.text = element_text(colour = "black",face = "bold",size=15),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=15),legend.title = element_text(colour = "black",face = "bold",size=15) ,strip.text = element_text(colour = "black",face = "bold",size=15))


svy.WHOrisk <- dat %>% 
            filter(is.na(WHOrisk)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk,
                        variables = c(WHOrisk,
                                      Country, countryRegionclass, region,
                                      sex, age_10yr,
                                      bmicat))
dat_avWHOrisk <- svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% summarise(avrisk=survey_quantile(WHOrisk,0.5,vartype="ci")) %>% mutate(medci=paste(sprintf("%.1f",avrisk_q50)," (",sprintf("%.1f",avrisk_q50_low),", ",sprintf("%.1f",avrisk_q50_upp),")",sep=""))%>% select(-avrisk_q50,-avrisk_q50_low,-avrisk_q50_upp)
dat_avWHOrisk<-dat_avWHOrisk %>% arrange(desc(region), factor(Country, levels = levels(orderCountry))) %>% mutate(sex=factor(sex,labels = c("Men","Women")))

p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_1", "meanWHOrisk_2", "meanWHOrisk_3", "meanWHOrisk_4", "meanWHOrisk_5"),labels=c("<5%","5 - <10%","10 - <20%","20 - <30%",">=30%"))))) + geom_bar(stat="identity",width = 0.9)+labs(y="Percent",x="")+geom_text(data=dat_avWHOrisk,aes(x=f(Country),y=14,label=medci,fontface=2),inherit.aes=FALSE,size=4)+coord_flip()+guides(fill=guide_legend(title="10-year \nCVD Risk"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=10))+ scale_fill_manual(values=c("#990000","orange","yellow","darkgreen","green"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=15),axis.text = element_text(colour = "black",face = "bold",size=15),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=15),legend.title = element_text(colour = "black",face = "bold",size=15) ,strip.text = element_text(colour = "black",face = "bold",size=15))

################  age-standardize

dat_agestan<-dat[dat$age_5yr!=11,]#drop 65-74 agegrp
numwomenglobal_GBD2017<-with(popdat,popdat[location_name=="Global" & sex_name =="Female" & age_group_name %in% c("30 to 34","35 to 39","40 to 44","45 to 49","50 to 54", "55 to 59","60 to 64") & year_id==2017,]$val)
wF_GBD2017<-numwomenglobal_GBD2017/sum(numwomenglobal_GBD2017)
nummenglobal_GBD2017<-with(popdat,popdat[location_name=="Global" & sex_name =="Male" & age_group_name %in% c("30 to 34","35 to 39","40 to 44","45 to 49","50 to 54", "55 to 59","60 to 64") & year_id==2017,]$val)
wM_GBD2017<-nummenglobal_GBD2017/sum(nummenglobal_GBD2017)
agestan<-rbind(data.frame(w_GBD2017=wF_GBD2017,sex=1,age_5yr=4:10),data.frame(w_GBD2017=wM_GBD2017,sex=0,age_5yr=4:10))

dat_agestan<-left_join(dat_agestan,agestan,by=c("sex","age_5yr"))
dat_agestan$weq_CVDrisk_agestan<-dat_agestan$weq_CVDrisk*dat_agestan$w_GBD2017
#dat_agestan<-dat_agestan %>% group_by(Country) %>% mutate(sumwt=sum(weq_CVDrisk_agestan))
#dat_agestan$weq_CVDrisk_agestan<-dat_agestan$weq_CVDrisk_agestan/dat_agestan$sumwt


svy.WHOrisk <- dat_agestan %>% 
            filter(is.na(WHOrisk)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk_agestan,
                        variables = c(WHOrisk,WHOrisk5cat,
                                      Country, countryRegionclass, region,
                                      sex, age_10yr,
                                      bmicat))
dat_WHOrisk <- svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% 
  summarize(meanWHOrisk_1 = survey_mean(WHOrisk5cat=="<5", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(WHOrisk5cat=="5-<10", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_3 = survey_mean(WHOrisk5cat=="10-<20", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_4 = survey_mean(WHOrisk5cat=="20-<30", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_5 = survey_mean(WHOrisk5cat==">=30", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se","meanWHOrisk_3_se","meanWHOrisk_4_se","meanWHOrisk_5_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_5),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_5)

f <- function(x) factor(x, levels = unique(x))
orderCountry<-dat_WHOrisk[dat_WHOrisk$RiskCategory=="meanWHOrisk_1" & dat_WHOrisk$sex==0,] %>% arrange(desc(region),desc(Pct)) %>% select(Country) %>%  .$Country %>% f(.)
dat_WHOrisk<-dat_WHOrisk %>%
   arrange(desc(region), factor(Country, levels = levels(orderCountry)))


p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_5", "meanWHOrisk_4", "meanWHOrisk_3", "meanWHOrisk_2", "meanWHOrisk_1"),labels=c(">=30%","20 - <30%","10 - <20%","5 - <10%","<5%"))))) + geom_bar(stat="identity",width = 0.9)+labs(y="Percent",x="")+coord_flip()+guides(fill=guide_legend(title="10-year \nCVD Risk"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=10))+ scale_fill_manual(values=c("green","darkgreen","yellow","orange","#990000"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=15),axis.text = element_text(colour = "black",face = "bold",size=15),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=15),legend.title = element_text(colour = "black",face = "bold",size=15) ,strip.text = element_text(colour = "black",face = "bold",size=15))


svy.WHOrisk <- dat_agestan %>% 
            filter(is.na(WHOrisk)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk_agestan,
                        variables = c(WHOrisk,WHOrisk5cat,
                                      Country, countryRegionclass, region,
                                      sex, age_10yr,
                                      bmicat))
dat_avWHOrisk <- svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% summarise(avrisk=survey_quantile(WHOrisk,0.5,vartype="ci"),meanWHOrisk_1 = survey_mean(WHOrisk5cat=="<5", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(WHOrisk5cat=="5-<10", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_3 = survey_mean(WHOrisk5cat=="10-<20", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_4 = survey_mean(WHOrisk5cat=="20-<30", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_5 = survey_mean(WHOrisk5cat==">=30", proportion=TRUE, prop_method="asin"))

################### text in manuscript   #####################
avCVDrisk_men<-round(c(quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.5),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.25),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.75),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0])),1)

avCVDrisk_women<-round(c(quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.5),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.25),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.75),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1])),1)

#low (<5%), low (5 -<10%), medium (10 - <20%), high (20- <30%) and very high risk (???30%)
#male
round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0])),1)

#female
round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1])),1)



dat_avWHOrisk <- dat_avWHOrisk %>% mutate(medci=paste(sprintf("%.1f",avrisk_q50)," (",sprintf("%.1f",avrisk_q50_low),", ",sprintf("%.1f",avrisk_q50_upp),")",sep=""))%>% select(-avrisk_q50,-avrisk_q50_low,-avrisk_q50_upp)
dat_avWHOrisk<-dat_avWHOrisk %>% arrange(desc(region), factor(Country, levels = levels(orderCountry))) %>% mutate(sex=factor(sex,labels = c("Men","Women")))

levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="India"]<-"India*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Ecuador"]<-"Ecuador*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Namibia"]<-"Namibia*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Albania"]<-"Albania*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Kyrgyz Republic"]<-"Kyrgyz Republic*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Lesotho"]<-"Lesotho*"

levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="India"]<-"India*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Ecuador"]<-"Ecuador*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Namibia"]<-"Namibia*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Albania"]<-"Albania*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Kyrgyz Republic"]<-"Kyrgyz Republic*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Lesotho"]<-"Lesotho*"

p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_5", "meanWHOrisk_4", "meanWHOrisk_3", "meanWHOrisk_2", "meanWHOrisk_1"),labels=c(paste("\u2265","30%",sep=""),"20 to <30%","10 to <20%","5 to <10%","<5%"))))) + geom_bar(stat="identity",width = 0.8)+labs(y="Percent",x="")+geom_text(data=dat_avWHOrisk,aes(x=f(Country),y=90,label=medci,fontface=2),inherit.aes=FALSE,size=3)+coord_flip()+guides(fill=guide_legend(title="10-year \nCVD Risk"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=12))+ scale_fill_manual(values=c("darkolivegreen3","yellow","orange","red","#990000"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=12),axis.text = element_text(colour = "black",face = "bold",size=10),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=12),legend.title = element_text(colour = "black",face = "bold",size=12) ,strip.text = element_text(colour = "black",face = "bold",size=12))
```

```{r : figure 1 for 40-plus}
dat_agestan<-dat[dat$age_5yr!=11 & dat$age_5yr!=4 & dat$age_5yr!=5,]#drop 65-74, 30-34 and 35-39 agegrp
numwomenglobal_GBD2017<-with(popdat,popdat[location_name=="Global" & sex_name =="Female" & age_group_name %in% c("40 to 44","45 to 49","50 to 54", "55 to 59","60 to 64") & year_id==2017,]$val)
wF_GBD2017<-numwomenglobal_GBD2017/sum(numwomenglobal_GBD2017)
nummenglobal_GBD2017<-with(popdat,popdat[location_name=="Global" & sex_name =="Male" & age_group_name %in% c("40 to 44","45 to 49","50 to 54", "55 to 59","60 to 64") & year_id==2017,]$val)
wM_GBD2017<-nummenglobal_GBD2017/sum(nummenglobal_GBD2017)
agestan<-rbind(data.frame(w_GBD2017=wF_GBD2017,sex=1,age_5yr=6:10),data.frame(w_GBD2017=wM_GBD2017,sex=0,age_5yr=6:10))

dat_agestan<-left_join(dat_agestan,agestan,by=c("sex","age_5yr"))
dat_agestan$weq_CVDrisk_agestan<-dat_agestan$weq_CVDrisk*dat_agestan$w_GBD2017
#dat_agestan<-dat_agestan %>% group_by(Country) %>% mutate(sumwt=sum(weq_CVDrisk_agestan))
#dat_agestan$weq_CVDrisk_agestan<-dat_agestan$weq_CVDrisk_agestan/dat_agestan$sumwt


svy.WHOrisk <- dat_agestan %>% 
            filter(is.na(WHOrisk)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk_agestan,
                        variables = c(WHOrisk,WHOrisk5cat,
                                      Country, countryRegionclass, region,
                                      sex, age_10yr,
                                      bmicat))
dat_WHOrisk <- svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% 
  summarize(meanWHOrisk_1 = survey_mean(WHOrisk5cat=="<5", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(WHOrisk5cat=="5-<10", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_3 = survey_mean(WHOrisk5cat=="10-<20", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_4 = survey_mean(WHOrisk5cat=="20-<30", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_5 = survey_mean(WHOrisk5cat==">=30", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se","meanWHOrisk_3_se","meanWHOrisk_4_se","meanWHOrisk_5_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_5),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_5)

f <- function(x) factor(x, levels = unique(x))
orderCountry<-dat_WHOrisk[dat_WHOrisk$RiskCategory=="meanWHOrisk_1" & dat_WHOrisk$sex==0,] %>% arrange(desc(region),desc(Pct)) %>% select(Country) %>%  .$Country %>% f(.)
dat_WHOrisk<-dat_WHOrisk %>%
   arrange(desc(region), factor(Country, levels = levels(orderCountry)))


p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_5", "meanWHOrisk_4", "meanWHOrisk_3", "meanWHOrisk_2", "meanWHOrisk_1"),labels=c(">=30%","20 - <30%","10 - <20%","5 - <10%","<5%"))))) + geom_bar(stat="identity",width = 0.9)+labs(y="Percent",x="")+coord_flip()+guides(fill=guide_legend(title="10-year \nCVD Risk"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=10))+ scale_fill_manual(values=c("green","darkgreen","yellow","orange","#990000"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=15),axis.text = element_text(colour = "black",face = "bold",size=15),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=15),legend.title = element_text(colour = "black",face = "bold",size=15) ,strip.text = element_text(colour = "black",face = "bold",size=15))


svy.WHOrisk <- dat_agestan %>% 
            filter(is.na(WHOrisk)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk_agestan,
                        variables = c(WHOrisk,WHOrisk5cat,
                                      Country, countryRegionclass, region,
                                      sex, age_10yr,
                                      bmicat))
dat_avWHOrisk <- svy.WHOrisk %>%
  group_by(region,Country,sex)  %>% summarise(avrisk=survey_quantile(WHOrisk,0.5,vartype="ci"),meanWHOrisk_1 = survey_mean(WHOrisk5cat=="<5", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(WHOrisk5cat=="5-<10", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_3 = survey_mean(WHOrisk5cat=="10-<20", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_4 = survey_mean(WHOrisk5cat=="20-<30", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_5 = survey_mean(WHOrisk5cat==">=30", proportion=TRUE, prop_method="asin"))

################### text in manuscript   #####################
avCVDrisk_men<-round(c(quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.5),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.25),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.75),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==0])),1)

avCVDrisk_women<-round(c(quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.5),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.25),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.75),quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$avrisk_q50[dat_avWHOrisk$sex==1])),1)

#low (<5%), low (5 -<10%), medium (10 - <20%), high (20- <30%) and very high risk (???30%)
#male
round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==0])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0]),quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0],0.25),range(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==0])),1)

#female
round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_1[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_2[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_3[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_4[dat_avWHOrisk$sex==1])),1)

round(100*c(quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1]),quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1],0.75)-quantile(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1],0.25),range(dat_avWHOrisk$meanWHOrisk_5[dat_avWHOrisk$sex==1])),1)



dat_avWHOrisk <- dat_avWHOrisk %>% mutate(medci=paste(sprintf("%.1f",avrisk_q50)," (",sprintf("%.1f",avrisk_q50_low),", ",sprintf("%.1f",avrisk_q50_upp),")",sep=""))%>% select(-avrisk_q50,-avrisk_q50_low,-avrisk_q50_upp)
dat_avWHOrisk<-dat_avWHOrisk %>% arrange(desc(region), factor(Country, levels = levels(orderCountry))) %>% mutate(sex=factor(sex,labels = c("Men","Women")))

levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="India"]<-"India*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Ecuador"]<-"Ecuador*"
#levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Namibia"]<-"Namibia*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Albania"]<-"Albania*"
#levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Kyrgyz Republic"]<-"Kyrgyz Republic*"
levels(dat_WHOrisk$Country)[levels(dat_WHOrisk$Country)=="Lesotho"]<-"Lesotho*"

levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="India"]<-"India*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Ecuador"]<-"Ecuador*"
#levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Namibia"]<-"Namibia*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Albania"]<-"Albania*"
#levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Kyrgyz Republic"]<-"Kyrgyz Republic*"
levels(dat_avWHOrisk$Country)[levels(dat_avWHOrisk$Country)=="Lesotho"]<-"Lesotho*"

p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_5", "meanWHOrisk_4", "meanWHOrisk_3", "meanWHOrisk_2", "meanWHOrisk_1"),labels=c(paste("\u2265","30%",sep=""),"20 to <30%","10 to <20%","5 to <10%","<5%"))))) + geom_bar(stat="identity",width = 0.8)+labs(y="Percent",x="")+geom_text(data=dat_avWHOrisk,aes(x=f(Country),y=90,label=medci,fontface=2),inherit.aes=FALSE,size=3)+coord_flip()+guides(fill=guide_legend(title="10-year \nCVD Risk"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=12))+ scale_fill_manual(values=c("darkolivegreen3","yellow","orange","red","#990000"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=12),axis.text = element_text(colour = "black",face = "bold",size=10),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=12),legend.title = element_text(colour = "black",face = "bold",size=12) ,strip.text = element_text(colour = "black",face = "bold",size=12))
```

```{r : Figure 2-multivariable regression CVD risk (medium-high) ~ edu+ WQ + marital status + empl}
dat <- dat %>% 
  mutate(poor = ifelse(wealth_quintile ==1|wealth_quintile ==2, 1, 0))
dat <- dat %>% 
  mutate(working = ifelse(working ==0, 0, ifelse(working==1,1,NA)))

md<-dat %>% group_by(Country) %>% summarise_at(.vars= vars( edulte4yr,poor,wpartner,working), funs(miss=sum(is.na(.))))%>% left_join(.,dat%>% group_by(Country) %>% summarise(n=n(),n.complete=sum(!is.na(edulte4yr) & !is.na(poor) & !is.na(wpartner) & !is.na(working)  ))) 
#%>% mutate_at(.vars= vars( edulte4yr_miss, poor_miss,wpartner_miss,working_miss,n.complete),funs(paste(.," (",sprintf("%.1f",.*100/n),")",sep=""))) 
countrylist<-dat %>% group_by(region,Country) %>% summarise(n=n()) %>% pull(Country)
est<-data.frame(country=rep(countrylist,each=4),rf=rep(c("edulte4yr", "poor","wpartner","working"), length(countrylist)),setNames(data.frame(matrix(NA,nrow=4*length(countrylist),ncol=2)),c("beta","se")))
countrylist_rev<-countrylist[!(countrylist %in% c("Uganda",unique(with(md[md$n.complete==0,],paste(Country)))))]
i<-0
for (country in countrylist_rev){
  i<-grep(country,countrylist)
  print(i)
  anadat<-na.omit(dat[dat$Country == country,c("Country","WHOrisk","edulte4yr","poor","female_lab","psu_num","weq_CVDrisk","wpartner","working")])
  if (dim(anadat)[1]>0){
    if (!(country %in% c("Uganda","Grenada","Morocco"))){
      if (sum(table(anadat$female_lab)<5)==0 & sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0) {m<-lmer(log(WHOrisk)~female_lab+I(edulte4yr==0)+I(poor==0)+wpartner+working+(1 | psu_num),data = anadat, weights = weq_CVDrisk)
est[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[-(1:2),1:2]
if (isSingular(m, tol = 1e-05)) print(paste ("singular fit - F&M - ",country)) 
    } 
        
  }
  else {
    if (sum(table(anadat$female_lab)<5)==0 & sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0) {m<-lm(log(WHOrisk)~female_lab+I(edulte4yr==0)+I(poor==0)+wpartner+working,data = anadat, weights = weq_CVDrisk)
est[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[-(1:2),1:2]
    } 
  }
}}
estval<-data.frame(est,pctchange=round(100*(exp(est$beta)-1),1),pctchange_l=round(100*(exp(est$beta-1.96*est$se)-1),1),pctchange_u=round(100*(exp(est$beta+1.96*est$se)-1),1))

####   plot 
regcoef <- read.xlsx("TabsFigs.xlsx",sheet="multi_all",colNames=TRUE)
dim(regcoef)
regcoef<-regcoef[!(regcoef$country %in% unique(regcoef$country[ regcoef$country%in% dat$Country & is.na(regcoef$beta_F)])),]
dim(regcoef)

fpdat<-regcoef[regcoef$level=="edulte4yr" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
xtcks_edu<-c(-90,0,60,120)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,8),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks_edu,align = c("l", "r"),mar = unit(c(-1,1,-3,0), "mm"),xlim=c(-90,120))

##  middle/richer/richest
fpdat<-regcoef[regcoef$level=="poor" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
#xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,8),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks_edu,align = c("l", "r"),mar = unit(c(-1,1,-3,0), "mm"),xlim=c(-90,120))

##  married
fpdat<-regcoef[regcoef$level=="wpartner" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,8),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks_edu,align = c("l", "r"),mar = unit(c(-1,1,-3,0), "mm"),xlim=c(-90,120))

# working
fpdat<-regcoef[regcoef$level=="working" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
xtcks<-c(-66,0,9)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,8),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks_edu,align = c("l", "r"),mar = unit(c(-1,1,-3,0), "mm"),xlim=c(-90,120))
```

```{r : figure 3 - % taking bp meds among ones indicated for meds}
#	10-yr CVD risk>=30%  OR  (SBP >=160 mm Hg  OR  DBP >=100 mm Hg)  
#	20% >=10-yr CVD risk < 30%    AND  (SBP >=140 mm Hg  OR  DBP >=90 mm Hg)
# Others

#guideline indicated - (10-yr CVD risk >=30%)  OR  (>=160 mm Hg  OR  DBP >=100 mm Hg)  OR ((20% >=10-yr CVD risk < 30%)    AND  (SBP >=140 mm Hg  OR  DBP >=90 mm Hg))

#########################################################################################
#number in response to reviewers' comment: A key reason for including people aged 30-39 years is that overall xx% of people taking BP medication are in this age range
dat_1<-dat
dat_1$age30to39<-with(dat_1,ifelse(age<40,1,0))
trt_svy <- dat_1 %>% filter(is.na(bpms_hypt_med)==F &  !is.na(age30to39)) %>% 
            as_survey_design(stratum = stratum_num,ids = psu_num, 
                        weights = weq_CVDrisk,variables = c(region,female_lab,Country,age30to39,bpms_hypt_med))
table_trt0 <- trt_svy%>% subset(bpms_hypt_med==1) %>% 
  group_by(region,Country,female_lab,bpms_hypt_med) %>% 
  summarize(trtsvy = survey_mean(age30to39,proportion=TRUE, prop_method="asin", 
                                vartype = "ci"))  %>% mutate(trtsvy=100*trtsvy,trtsvy_low=100*trtsvy_low,trtsvy_upp=100*trtsvy_upp) 
quantile(table_trt0$trtsvy)
tapply(table_trt0$trtsvy,table_trt0$female_lab,quantile,na.rm=T)
#########################################################################################


dat <- dat %>% 
  mutate(riskcat = ifelse(WHOrisk>=30 | sbp_avg_rev>=160 | dbp_avg_rev >=100, "elevated", ifelse(WHOrisk>=20 & (sbp_avg_rev>=140 | dbp_avg_rev >=90),"moderate","others")))
dat <- dat %>% 
  mutate(guideRx_0 = ifelse(riskcat=="elevated" | riskcat=="moderate", "indicated","NOTindicated"))
dat <- dat %>% 
  mutate(guideRx = ifelse(WHOrisk>=30 | sbp_adjtrt>=160 | dbp_adjtrt >=100 | ((WHOrisk>=20 & WHOrisk<30) & (sbp_adjtrt>=140 | dbp_adjtrt >=90)), "indicated","NOTindicated"))

dat$guideRx_lab<-factor(dat$guideRx,levels=c("NOTindicated","indicated"))

trt_svy <- dat %>% filter(is.na(bpms_hypt_med)==F & is.na(guideRx)==F ) %>% 
            as_survey_design(stratum = stratum_num,ids = psu_num, 
                        weights = weq_CVDrisk,variables = c(region,female_lab,Country,guideRx, 
                                    bpms_hypt_med))
table_trt0 <-trt_svy%>% 
  group_by(region,Country,female_lab,guideRx) %>% 
  summarize(trtsvy = survey_mean(bpms_hypt_med,proportion=TRUE, prop_method="asin", 
                                vartype = "ci"))  %>% mutate(trtsvy=100*trtsvy,trtsvy_low=100*trtsvy_low,trtsvy_upp=100*trtsvy_upp) 

ctl_svy <- dat %>% filter( is.na(bpms_hypt_med_ctl)==F & is.na(guideRx)==F ) %>% 
            as_survey_design(stratum = stratum_num,ids = psu_num, 
                        weights = weq_CVDrisk,variables = c(region,female_lab,Country,guideRx, 
                                      bpms_hypt_med_ctl))
table_ctl0 <-ctl_svy%>% 
  group_by(region,Country,female_lab,guideRx) %>% 
  summarize(ctlsvy = survey_mean(bpms_hypt_med_ctl,proportion=TRUE, prop_method="asin", 
                                vartype = "ci"))  %>% mutate(ctlsvy=100*ctlsvy,ctlsvy_low=100*ctlsvy_low,ctlsvy_upp=100*ctlsvy_upp)
#table_ctl0<-left_join(table_ctl0,unique(dat[,c("Country","region1")]))

table_trt1<-table_trt0 %>% mutate(val=trtsvy,lower=trtsvy_low,upper=trtsvy_upp,which="trt") %>% select(-one_of(c("trtsvy","trtsvy_low","trtsvy_upp")))
table_ctl1<-table_ctl0 %>% mutate(val=ctlsvy,lower=ctlsvy_low,upper=ctlsvy_upp,which="ctl") %>% select(-one_of(c("ctlsvy","ctlsvy_low","ctlsvy_upp")))
table_trtctl1<-bind_rows(table_trt1,table_ctl1)

#Treated
plotdat<-table_trtctl1 %>% filter(which=="trt" & guideRx=="indicated")
f <- function(x) factor(x, levels = unique(x))
orderCountry<-plotdat[plotdat$guideRx=="indicated" & plotdat$female_lab=="Male",] %>% arrange(desc(region),desc(val)) %>% select(Country) %>%  .$Country %>% f(.)
plotdat<-plotdat %>% arrange(desc(region), factor(Country, levels = levels(orderCountry)))

#### ### ### ###   numbers in manuscript   #### ### ### ###
BPmed_indicatedMen<- plotdat$val[plotdat$guideRx=="indicated" & plotdat$female_lab=="Male"]
avBPmed_indicatedMen<-round(c(quantile(BPmed_indicatedMen),quantile(BPmed_indicatedMen,0.75)-quantile(BPmed_indicatedMen,0.25),range(BPmed_indicatedMen,na.rm=T)),1)

BPmed_indicatedWomen<- plotdat$val[plotdat$guideRx=="indicated" & plotdat$female_lab=="Female"]
avBPmed_indicatedWomen<-round(c(quantile(BPmed_indicatedWomen),quantile(BPmed_indicatedWomen,0.75)-quantile(BPmed_indicatedWomen,0.25),range(BPmed_indicatedWomen,na.rm=T)),1)
#### ### ### ### #### ### ### ### #### ### ### ### #### ### ### ###

p<-plotdat %>% ggplot(aes(x=f(Country),y=val))+geom_bar(stat="identity", position=position_dodge())+geom_errorbar(aes(ymin=lower, ymax=upper), size=.3,    width=.2,       position=position_dodge(.9)) +  labs(x="", y="")+ theme(panel.background = element_rect(fill = NA),axis.text = element_text(colour = "black",face = "bold",size=12),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, 0, "cm"),strip.background.x = element_blank(),strip.background.y = element_rect(fill = "grey95"),legend.position = "bottom",legend.margin=margin(0,0,0,0),legend.box.margin=margin(-22,-1,1,-1),strip.text = element_text(size = 11,face = "bold"),legend.text=element_text(size=12),legend.title=element_text(size=12) )+coord_flip()+ guides(color = FALSE,fill = guide_legend(reverse = TRUE))+facet_grid( region~fct_rev(factor(female_lab,labels = c("Women","Men"))),scales="free", space="free",labeller = label_wrap_gen(width=20))

##plot with summary
plotdat_sum<-plotdat %>% group_by(region,female_lab) %>% summarise(country=unique(region),val_med=quantile(val,0.5),val_lower=quantile(val,0.25),val_upper=quantile(val,0.75),which="trt")
plotdat<-read.xlsx("TabsFigs.xlsx",sheet="BPmed_overall",colNames=TRUE)
p<-plotdat %>% ggplot(aes(x=f(Country),y=val))+geom_bar(stat="identity", position=position_dodge())+geom_errorbar(aes(ymin=lower, ymax=upper), size=.3,    width=.2,       position=position_dodge(.9)) +  labs(x="", y="")+ theme(panel.background = element_rect(fill = NA),axis.text = element_text(colour = "black",face = "bold",size=12),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, 0, "cm"),strip.background.x = element_blank(),strip.background.y = element_rect(fill = "grey95"),legend.position = "bottom",legend.margin=margin(0,0,0,0),legend.box.margin=margin(-22,-1,1,-1),strip.text = element_text(size = 11,face = "bold"),legend.text=element_text(size=12),legend.title=element_text(size=12) )+coord_flip()+ guides(color = FALSE,fill = guide_legend(reverse = TRUE))+facet_grid( region~fct_rev(factor(female_lab,labels = c("Women","Men"))),scales="free", space="free",labeller = label_wrap_gen(width=20))
p  

#plot with RRs

## RRs
countrylist<-dat %>% group_by(region,Country) %>% summarise(n=n()) %>% pull(Country)
#countrylist<-unique(dat$Country)[(unique(dat$Country) %in% newdata$Country)]

est_F<-est_M<-data.frame(country=rep(countrylist,each=2),edulevel=rep(c("NOTindicated","guideRx_labindicated"),length(countrylist)),setNames(data.frame(matrix(NA,nrow=2*length(countrylist),ncol=4)),c("n","beta","se","pval")))
i<-0
for (country in countrylist){
  i<-i+1
  if (country !="Kyrgyz Republic") {
  anadat<-na.omit(dat[dat$Country == country,c("Country","guideRx_lab","bpms_hypt_med","female_lab","psu_num","weq_CVDrisk","pid")])
if (dim(anadat)[1]>0){
    est_F$n[(2*i-1):(2*i)]<-table(anadat$guideRx_lab[anadat$female_lab=="Female"])
    est_M$n[(2*i-1):(2*i)]<-table(anadat$guideRx_lab[anadat$female_lab=="Male"])
    num.zerocell_F<-with(anadat[anadat$female_lab=="Female",],sum(table(guideRx_lab,bpms_hypt_med)<5))
    num.zerocell_M<-with(anadat[anadat$female_lab=="Male",],sum(table(guideRx_lab,bpms_hypt_med)<5))
    
    if (num.zerocell_F>0)  print(paste("Female- ",country) )
    if (num.zerocell_M>0)  print(paste("Male- ",country) )
    if (num.zerocell_F==0) {m_F<-geeglm(formula = bpms_hypt_med~guideRx_lab,data=anadat[anadat$female_lab=="Female",], family  = poisson(link = "log"), id= pid, corstr  = "exchangeable",weights=weq_CVDrisk)
    est_F[(2*i),][match(rownames(summary(m_F)$coefficients)[-1],est_F[(2*i),2]),4:6]<-summary(m_F)$coefficients[-1,c(1:2,4)]}
    if (num.zerocell_M==0) {m_M<-geeglm(formula = bpms_hypt_med~guideRx_lab,data=anadat[anadat$female_lab=="Male",], family  = poisson(link = "log"), id= pid, corstr  = "exchangeable",weights=weq_CVDrisk)
   est_M[(2*i),][match(rownames(summary(m_M)$coefficients)[-1],est_M[(2*i),2]),4:6]<-summary(m_M)$coefficients[-1,c(1:2,4)]}}}
  
  if (country =="Kyrgyz Republic") {
  anadat<-na.omit(dat[dat$Country == country,c("Country","guideRx_lab","bpms_hypt_med","female_lab","psu_num","weq_CVDrisk","pid")])
if (dim(anadat)[1]>0){
    est_F$n[(2*i-1):(2*i)]<-table(anadat$guideRx_lab[anadat$female_lab=="Female"])
    num.zerocell_F<-with(anadat[anadat$female_lab=="Female",],sum(table(guideRx_lab,bpms_hypt_med)<5))
    
    if (num.zerocell_F>0)  print(paste("Female- ",country) )
    if (num.zerocell_F==0) {m_F<-geeglm(formula = bpms_hypt_med~guideRx_lab,data=anadat[anadat$female_lab=="Female",], family  = poisson(link = "log"), id      = pid, corstr  = "exchangeable",weights=weq_CVDrisk)
    est_F[(2*i),][match(rownames(summary(m_F)$coefficients)[-1],est_F[(2*i),2]),4:6]<-summary(m_F)$coefficients[-1,c(1:2,4)]}
    }}
   print(i)
  }
estpval_F<-data.frame(est_F,pctchange=round((exp(est_F$beta)),1),pctchange_l=round((exp(est_F$beta-1.96*est_F$se)),1),pctchange_u=round((exp(est_F$beta+1.96*est_F$se)),1))
estpval_M<-data.frame(est_M,pctchange=round((exp(est_M$beta)),1),pctchange_l=round((exp(est_M$beta-1.96*est_M$se)),1),pctchange_u=round((exp(est_M$beta+1.96*est_M$se)),1))

levels(estpval_F$edulevel)[levels(estpval_F$edulevel)=="guideRx_labindicated"]<-"indicated"
levels(estpval_M$edulevel)[levels(estpval_M$edulevel)=="guideRx_labindicated"]<-"indicated"

# text in manuscript
#univariate relative risk (RR (95% CI)) of receiving BP medication ranged from 1.6 (x.x - x.x) in Costa Rica to 9.9 (xx.x - xx.x) in Burkina Faso
estpval_F[which.max(estpval_F$pctchange),]
estpval_F[which.min(estpval_F$pctchange),]
#For men, the corresponding RRs ranged from 2.3 (x.x - x.x) in Granada to 11.1 (x.x - xx.x) in Tanzania. 
estpval_M[which.max(estpval_M$pctchange),]
estpval_M[which.min(estpval_M$pctchange),]

regcoef_edu <- read.xlsx("TabsFigs.xlsx",sheet="BPmed_rev",colNames=TRUE)

names(regcoef_edu)[1:2]<-c("Country", "guideRx")
plotdatF<-left_join(plotdat[plotdat$female_lab=="Female",],regcoef_edu[,c("Country", "guideRx","pctchange_F","pctchange_l_F", "pctchange_u_F")])
plotdatF <-plotdatF %>% mutate(RR=ifelse(is.na(pctchange_F),NA,paste(sprintf("%.1f",pctchange_F)," (", sprintf("%.1f",pctchange_l_F),", ", sprintf("%.1f",pctchange_u_F),")",sep=""))) %>%  select(region,Country,female_lab,guideRx,val,lower,upper,which,RR)
#names(plotdatF)[dim(plotdatF)[2]]<-"RR"
plotdatM<-left_join(plotdat[plotdat$female_lab=="Male",],regcoef_edu[,c("Country", "guideRx","pctchange_M","pctchange_l_M", "pctchange_u_M")])
plotdatM <-plotdatM %>% mutate(RR=ifelse(is.na(pctchange_M),NA,paste(sprintf("%.1f",pctchange_M)," (", sprintf("%.1f",pctchange_l_M),", ", sprintf("%.1f",pctchange_u_M),")",sep=""))) %>%  select(region,Country,female_lab,guideRx,val,lower,upper,which,RR)
#names(plotdatM)[dim(plotdatM)[2]]<-"RR"
plotdat1<-rbind(plotdatF,plotdatM)
#plotdat1$RR<-with(plotdat1,ifelse(is.na(RR),NA,sprintf("%.1f",RR)))
#plotdat1$RR_hgh<-with(plotdat1,ifelse(is.na(RR),NA,ifelse(WHOriskcat=="High",sprintf("%.1f",RR),NA)))

p<-plotdat1 %>% ggplot(aes(x=f(Country),y=val))+geom_bar(stat="identity", position=position_dodge())+geom_text(aes(y=107,label=RR ), size = 3.5) +  scale_y_continuous(limits = c(0,110),breaks = c(0,10,20,30,40,50,60,70,80, 90, 100)) +geom_errorbar(aes(ymin=lower, ymax=upper), size=.3,    width=.2,       position=position_dodge(.9)) +  labs(y="Percent",x="")+ theme(panel.background = element_rect(fill = NA),axis.text.x = element_text(colour = "black",face = "bold",size=7),axis.text.y = element_text(colour = "black",face = "bold",size=10),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, 0, "cm"),strip.background.x = element_blank(),strip.background.y = element_rect(fill = "grey95"),legend.position = "bottom",legend.margin=margin(0,0,0,0),legend.box.margin=margin(-22,-1,1,-1),strip.text = element_text(size = 11,face = "bold"),legend.text=element_text(size=12),legend.title=element_text(size=12) )+coord_flip()+ guides(color = FALSE,fill = guide_legend(reverse = TRUE))+facet_grid( region~fct_rev(factor(female_lab,labels = c("Women","Men"))),scales="free", space="free",labeller = label_wrap_gen(width=12))
p<-p + theme(  panel.background = element_rect(fill = "white",
                                colour = "white",
                                size = 0.5, linetype = "solid"),
  panel.grid.major.x = element_line(size = 0.5, linetype = 'dashed',colour = "gray"), 
  panel.grid.major.y = element_blank()
  #panel.grid.minor = element_line(size = 0.25, linetype = 'solid',   colour = "gray")
  )
p

#cleaveland plot - Pmed reviewer

p<-plotdat1 %>% ggplot(aes(x=f(Country),y=val))+geom_point(size=1)+geom_text(aes(y=90,label=RR ), size = 3,hjust = 0) +  scale_y_continuous(limits = c(0,137),breaks = c(0,10,20,30,40,50,60,70,80,90)) +geom_errorbar(aes(ymin=lower, ymax=upper), size=.3,    width=.2,       position=position_dodge(.9)) +  labs(y="Percent",x="")+ theme(panel.background = element_rect(fill = NA),axis.text.x = element_text(colour = "black",face = "bold",size=7),axis.text.y = element_text(colour = "black",face = "bold",size=10),panel.border = element_rect(fill = NA),plot.margin = margin(0,0.2 , 0, 0, "cm"),strip.background.x = element_blank(),strip.background.y = element_rect(fill = "grey95"),legend.position = "bottom",legend.margin=margin(0,0,0,0),legend.box.margin=margin(-22,-1,1,-1),strip.text = element_text(size = 11,face = "bold"),legend.text=element_text(size=12),legend.title=element_text(size=12) )+coord_flip()+ guides(color = FALSE,fill = guide_legend(reverse = TRUE))+facet_grid( region~fct_rev(factor(female_lab,labels = c("Women","Men"))),scales="free", space="free",labeller = label_wrap_gen(width=12))
p<-p + theme(  panel.background = element_rect(fill = "white",
                                colour = "white",
                                size = 0.5, linetype = "solid"),
  panel.grid.major.x = element_line(size = 0.5, linetype = 'dashed',colour = "gray"), 
  panel.grid.major.y = element_blank()
  #panel.grid.minor = element_line(size = 0.25, linetype = 'solid',   colour = "gray")
  )
p


```

```{r : Figure 4: elevated risk distribution among the ones taking BP meds}
# elevated risk acc to WHO/ISH guidelines - high CVD risk and/or high bp
#	High risk and indicated for BP Rx (10-yr CVD risk ??? 30%  OR  SBP ???160 mm Hg  OR  DBP ???110 mm Hg)
#	Moderate risk and indicated for BP Rx (20% ???10-yr CVD risk < 30% AND (SBP???140 mm Hg  OR  DBP ???90 mm Hg))
#	All others (low to moderate risk and not indicated for Rx)

svy.WHOrisk <- dat %>% 
            filter(is.na(guideRx)==F) %>% 
            as_survey_design(stratum = stratum_num,
                        ids = c(psu_num), 
                        weights = weq_CVDrisk,
                        variables = c(WHOrisk,WHOriskcat, bpms_hypt_med,riskcat,guideRx_lab,
                                      Country, countryRegionclass, region,age, currsmok,bpms_hypt_med,sbp_avg_rev,bmi,bpms_hypt_med,
                                      sex, age_10yr,
                                      bmicat))
dat_WHOrisk <- svy.WHOrisk %>% filter(bpms_hypt_med==1)%>%
  group_by(region,Country,sex)  %>% 
  summarize(meanWHOrisk_1 = survey_mean(guideRx_lab=="NOTindicated", proportion=TRUE, prop_method="asin"),
            meanWHOrisk_2 = survey_mean(guideRx_lab=="indicated", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_2),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_2)

f <- function(x) factor(x, levels = unique(x))
orderCountry<-dat_WHOrisk[dat_WHOrisk$RiskCategory=="meanWHOrisk_1" & dat_WHOrisk$sex==0,] %>% arrange(desc(region),desc(Pct)) %>% select(Country) %>%  .$Country %>% f(.)
dat_WHOrisk<-dat_WHOrisk %>%
   arrange(desc(region), factor(Country, levels = levels(orderCountry)))

p<-dat_WHOrisk %>% mutate(sex=factor(sex,labels = c("Men","Women"))) %>% ggplot(aes(y = Pct, x = f(Country), fill = forcats::fct_rev(factor(RiskCategory,levels=c("meanWHOrisk_2", "meanWHOrisk_1"),labels=c("Yes","No"))))) + geom_bar(stat="identity",width = 0.9)+labs(y="Percent",x="")+coord_flip()+guides(fill=guide_legend(title="Indicated for\n medication"))+facet_grid(region~factor(sex),scales = "free", space = "free",labeller = label_wrap_gen(width=10))+scale_fill_manual(values=c("#990000","orange"))+theme(panel.background = element_rect(fill = NA),axis.title = element_text(colour = "black",face = "bold",size=15),axis.text = element_text(colour = "black",face = "bold",size=12),panel.border = element_rect(fill = NA),plot.margin = margin(-0.2,0.2 , 0, -0.45, "cm"),strip.background = element_blank(),legend.text =element_text(colour = "black",face = "bold",size=15),legend.title = element_text(colour = "black",face = "bold",size=15) ,strip.text = element_text(colour = "black",face = "bold",size=15))

# proportion of people taking BP medication who were not indicated for such based on their CVD risk status
dat_WHOrisk <- svy.WHOrisk %>% filter(bpms_hypt_med==1)%>%
    group_by(region,Country)  %>% 
    summarize(meanWHOrisk_1 = survey_mean(guideRx_lab=="NOTindicated", proportion=TRUE, prop_method="asin"),
              meanWHOrisk_2 = survey_mean(guideRx_lab=="indicated", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_2),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_2)

avNoIndicated<-round(c(quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"]),quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],0.75)-quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],0.25),range(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],na.rm=T)),1)
avNoIndicated

#men
dat_WHOrisk <- svy.WHOrisk %>% filter(bpms_hypt_med==1 & sex==0)%>%
    group_by(region,Country)  %>% 
    summarize(meanWHOrisk_1 = survey_mean(guideRx_lab=="NOTindicated", proportion=TRUE, prop_method="asin"),
              meanWHOrisk_2 = survey_mean(guideRx_lab=="indicated", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_2),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_2)

avNoIndicated_men<-round(c(quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"]),quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],0.75)-quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],0.25),range(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],na.rm=T)),1)


#women
dat_WHOrisk <- svy.WHOrisk %>% filter(bpms_hypt_med==1 & sex==1)%>%
    group_by(region,Country)  %>% 
    summarize(meanWHOrisk_1 = survey_mean(guideRx_lab=="NOTindicated", proportion=TRUE, prop_method="asin"),
              meanWHOrisk_2 = survey_mean(guideRx_lab=="indicated", proportion=TRUE, prop_method="asin")) %>% select(-one_of(c("meanWHOrisk_1_se","meanWHOrisk_2_se")))%>% mutate_at(vars(meanWHOrisk_1:meanWHOrisk_2),funs(.*100))%>% gather(.,key = RiskCategory, value = Pct, meanWHOrisk_1:meanWHOrisk_2)

avNoIndicated_women<-round(c(quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"]),quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],0.75)-quantile(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],0.25),range(dat_WHOrisk$Pct[dat_WHOrisk$RiskCategory=="meanWHOrisk_1"],na.rm=T)),1)

avNoIndicated_men
avNoIndicated_women

```

```{r : Figures S8 and S9 - SES and BP meds}
md<-dat %>% group_by(Country) %>% summarise_at(.vars= vars( edulte4yr,poor,wpartner,working,guideRx_lab,bpms_hypt_med), funs(miss=sum(is.na(.))))%>% left_join(.,dat%>%  group_by(Country) %>% summarise(n=n(),n.complete=sum(!is.na(edulte4yr) & !is.na(poor) & !is.na(wpartner) & !is.na(working)  & !is.na(guideRx_lab) &!is.na(bpms_hypt_med) ))) 
#%>% mutate_at(.vars= vars( edulte4yr_miss, poor_miss,wpartner_miss,working_miss,n.complete),funs(paste(.," (",sprintf("%.1f",.*100/n),")",sep=""))) 
countrylist<-dat %>% group_by(region,Country) %>% summarise(n=n()) %>% pull(Country)
est_NOTindicated<-data.frame(country=rep(countrylist,each=4),rf=rep(c("edulte4yr", "poor","wpartner","working"), length(countrylist)),setNames(data.frame(matrix(NA,nrow=4*length(countrylist),ncol=2)),c("beta","se")))
est_indicated<-data.frame(country=rep(countrylist,each=4),rf=rep(c("edulte4yr", "poor","wpartner","working"), length(countrylist)),setNames(data.frame(matrix(NA,nrow=4*length(countrylist),ncol=2)),c("beta","se")))

countrylist_rev<-countrylist[!(countrylist %in% c("Uganda","Tajikistan","Georgia",unique(with(md[md$n.complete==0,],paste(Country)))))]
i<-0
for (country in countrylist_rev){
  i<-grep(country,countrylist)
  print(i)
  anadat<-na.omit(dat[dat$Country == country,c("Country","guideRx_lab","bpms_hypt_med","edulte4yr","poor","female_lab","psu_num","weq_CVDrisk","wpartner","working","pid")])
  if (dim(anadat)[1]>0){
       if (sum(table(anadat$female_lab)<5)==0 & sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0 & sum(table(anadat$guideRx_lab )<5)==0) {
        #m<-geeglm(formula = bpms_hypt_med~female_lab+guideRx_lab+I(edulte4yr==0)+I(poor==0)+wpartner+working,data=anadat, family  = poisson(link = "log"), id= pid, corstr  = "exchangeable",weights=weq_CVDrisk)
        m<-geeglm(formula = bpms_hypt_med~female_lab+guideRx_lab+I(edulte4yr==0)+I(poor==0)+wpartner+working+guideRx_lab:I(edulte4yr==0)+guideRx_lab:I(poor==0)+guideRx_lab:wpartner+guideRx_lab:working,data=anadat, family  = poisson(link = "log"), id= pid, corstr  = "exchangeable",weights=weq_CVDrisk)
        est_NOTindicated[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[4:7,1:2]
        
        m<-geeglm(formula = bpms_hypt_med~female_lab+I(guideRx_lab=="NOTindicated")+I(edulte4yr==0)+I(poor==0)+wpartner+working+I(guideRx_lab=="NOTindicated"):I(edulte4yr==0)+I(guideRx_lab=="NOTindicated"):I(poor==0)+I(guideRx_lab=="NOTindicated"):wpartner+I(guideRx_lab=="NOTindicated"):working,data=anadat, family  = poisson(link = "log"), id= pid, corstr  = "exchangeable",weights=weq_CVDrisk)
        est_indicated[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[4:7,1:2]
#if (isSingular(m, tol = 1e-05)) print(paste ("singular fit - F&M - ",country)) 
    } 
  }}

estpval_NOTindicated<-data.frame(est_NOTindicated,pctchange=round((exp(est_NOTindicated$beta)),1),pctchange_l=round((exp(est_NOTindicated$beta-1.96*est_NOTindicated$se)),1),pctchange_u=round((exp(est_NOTindicated$beta+1.96*est_NOTindicated$se)),1))

estpval_indicated<-data.frame(est_indicated,pctchange=round((exp(est_indicated$beta)),1),pctchange_l=round((exp(est_indicated$beta-1.96*est_indicated$se)),1),pctchange_u=round((exp(est_indicated$beta+1.96*est_indicated$se)),1))

#########  plot - not indicated - S9
 
regcoef <- read.xlsx("TabsFigs.xlsx",sheet="medSES_NOTindicated",colNames=TRUE)
dim(regcoef)
regcoef<-regcoef[!(regcoef$country %in% unique(regcoef$country[ regcoef$country%in% dat$Country & is.na(regcoef$beta_F)])),]
dim(regcoef)

## prim school or above
fpdat<-regcoef[regcoef$level=="edulte4yr" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=c(0.1,11.3, 22.6),align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  middle/richer/richest
fpdat<-regcoef[regcoef$level=="poor" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=c(0.1 ,5.5 ,11.1),align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  married
fpdat<-regcoef[regcoef$level=="wpartner" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

# working
fpdat<-regcoef[regcoef$level=="working" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

#########  plot - indicated - Fig S8
 
regcoef <- read.xlsx("TabsFigs.xlsx",sheet="medSES_indicated",colNames=TRUE)
dim(regcoef)
regcoef<-regcoef[!(regcoef$country %in% unique(regcoef$country[ regcoef$country%in% dat$Country & is.na(regcoef$beta_F)])),]
dim(regcoef)

## prim school or above
fpdat<-regcoef[regcoef$level=="edulte4yr" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  middle/richer/richest
fpdat<-regcoef[regcoef$level=="poor" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  married
fpdat<-regcoef[regcoef$level=="wpartner" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

# working
fpdat<-regcoef[regcoef$level=="working" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],1,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]),1)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",sprintf("%.1f",pctchange_F)))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Risk Ratio       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,6),rep(TRUE,2),rep(FALSE,9),TRUE,rep(FALSE,14)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=1, cex=1, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=c(0.1,4.9, 9.8),align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))
```

```{r : table S4 missing data}
mdat<-dat1 %>% group_by(Country) %>% mutate(n_complete=sum(complete==1)) %>% filter(n_complete>0)
mdat<-filter (mdat,complete!=1 )
dim(mdat)

mdat$agegrp<-with(mdat,ifelse(age<35,"30-34",ifelse(age<45,"35-44",ifelse(age<55,"45-54",ifelse(age<65,"55-64","65+")))))
dat$agegrp<-with(dat,ifelse(age<35,"30-34",ifelse(age<45,"35-44",ifelse(age<55,"45-54",ifelse(age<65,"55-64","65+")))))

rbind(t(bind_rows(round(100*prop.table(table(mdat$female_lab,useNA="ifany")),1),
round(100*prop.table(table(dat$female_lab,useNA="ifany")),1))),
t(bind_rows(round(100*prop.table(table(mdat$agegrp,useNA="ifany")),1),
round(100*prop.table(table(dat$agegrp,useNA="ifany")),1))),
t(bind_rows(round(100*prop.table(table(mdat$bmicat_lab,useNA="ifany")),1),
round(100*prop.table(table(dat$bmicat_lab,useNA="ifany")),1))),
t(bind_rows(round(100*prop.table(table(mdat$currsmok,useNA="ifany")),1),
round(100*prop.table(table(dat$currsmok,useNA="ifany")),1))),
t(rbind(round(100*prop.table(table(mdat$edulte4yr,useNA="ifany")),1),
round(100*prop.table(table(dat$edulte4yr,useNA="ifany")),1))),
t(rbind(round(100*prop.table(table(mdat$wealth_quintile,useNA="ifany")),1),
round(100*prop.table(table(dat$wealth_quintile,useNA="ifany")),1))))

#####  weighted
# mdat weights taht sum up to 1 within mdat. Alternatively we can use data weights (weights_eq)
mdat$mweights<-mdat$w2
mdat<-mdat %>% group_by(Country) %>% mutate(num_obs_m=ifelse(!is.na(mweights),sum(mweights,na.rm=T),NA))
mdat$mweights_eq<-rep(NA,dim(mdat)[1])
mdat$mweights_eq=with(mdat,ifelse(!is.na(mweights),mweights/num_obs_m,NA))

mdat$stratum_num <- as.factor(mdat$stratum_num) 
mdat$psu_num <- as.factor(mdat$psu_num) 
mdat$hh_id <- as.factor(mdat$hh_id) 
mdat$weights_eq <- as.numeric(mdat$weights_eq) 
mdat$mweights_eq <- as.numeric(mdat$mweights_eq) 
mdat$Country <- as.factor(mdat$Country) 
# Set psu_num in those countries with only one PSU equal to hh_id
onlyonepsu <- c("China", "Grenada", "Kazakhstan", "Mexico", "Romania", "Seychelles", "Uganda")
mdat <- mdat %>% 
  mutate(psu_num = ifelse(Country %in% onlyonepsu, hh_id, psu_num))

tableS2.svy.1 <- mdat %>% 
            as_survey_design(stratum = stratum_num,
                        ids = psu_num, 
                        weights = mweights_eq,
                        variables = c(region,Country,sex,female_lab,age,currsmok,bpms_hypt_med,sbp_avg_rev,bmicat_lab,edulte4yr,wealth_quintile,agegrp))

tableS2.m.1<-tableS2.svy.1 %>% summarise(pctfem=survey_mean((female_lab=="Female" & !is.na(female_lab)),proportion=TRUE, prop_method="asin"),pctm=survey_mean((female_lab=="Male" & !is.na(female_lab)),proportion=TRUE, prop_method="asin"),pctsexmiss=survey_mean((is.na(female_lab)),proportion=TRUE, prop_method="asin"),pct30to34=survey_mean((agegrp=="30-34" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct35to44=survey_mean((agegrp=="35-44" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct45to54=survey_mean((agegrp=="45-54" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct55to64=survey_mean((agegrp=="55-64" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct65plus=survey_mean((agegrp=="65+" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"),pctagemiss=survey_mean(is.na(agegrp),proportion=TRUE, prop_method="asin"), pctunderwt=survey_mean((bmicat_lab=="<18.5" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"), pctnormal=survey_mean((bmicat_lab=="18.5-24.9" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"), pctoverwt=survey_mean((bmicat_lab=="25.0-29.9" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"),pctobese=survey_mean((bmicat_lab==">=30.0" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"),pctbmimiss=survey_mean(is.na(bmicat_lab),proportion=TRUE, prop_method="asin"),pctnonsmoker=survey_mean(currsmok==0 & !is.na(currsmok),proportion=TRUE, prop_method="asin"),pctsmoker=survey_mean(currsmok==1 & !is.na(currsmok),proportion=TRUE, prop_method="asin"),pctsmokemiss=survey_mean(is.na(currsmok),proportion=TRUE, prop_method="asin"),pctedultp=survey_mean(edulte4yr==1 & !is.na(edulte4yr),proportion=TRUE, prop_method="asin"),pcteduprimplus=survey_mean(edulte4yr==0 & !is.na(edulte4yr),proportion=TRUE, prop_method="asin"),pctedumiss=survey_mean(is.na(edulte4yr),proportion=TRUE, prop_method="asin"), pctWQ1=survey_mean(wealth_quintile==1 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQ2=survey_mean(wealth_quintile==2 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQ3=survey_mean(wealth_quintile==3 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"), pctWQ4=survey_mean(wealth_quintile==4 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQ5=survey_mean(wealth_quintile==5 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQmiss=survey_mean(is.na(wealth_quintile),proportion=TRUE, prop_method="asin")) %>% dplyr::select(pctfem,pctm,pctsexmiss,pct30to34,pct35to44,pct45to54,pct55to64,pct65plus,pctagemiss,pctunderwt,pctnormal,pctoverwt,pctobese,pctbmimiss,pctnonsmoker,pctsmoker,pctsmokemiss,pctedultp,pcteduprimplus,pctedumiss,pctWQ1,pctWQ2,pctWQ3,pctWQ4,pctWQ5,pctWQmiss) %>% mutate_all(funs( round(.* 100,1))) 
t(tableS2.m.1)


dat$stratum_num <- as.factor(dat$stratum_num) 
dat$psu_num <- as.factor(dat$psu_num) 
dat$hh_id <- as.factor(dat$hh_id) 
dat$weq_CVDrisk <- as.numeric(dat$weq_CVDrisk) 
dat$Country <- as.factor(dat$Country) 
# Set psu_num in those countries with only one PSU equal to hh_id
onlyonepsu <- c("China", "Grenada", "Kazakhstan", "Mexico", "Romania", "Seychelles", "Uganda")
dat <- dat %>% 
  mutate(psu_num = ifelse(Country %in% onlyonepsu, hh_id, psu_num))

tableS2.svy <- dat %>% 
            as_survey_design(stratum = stratum_num,
                        ids = psu_num, 
                        weights = weq_CVDrisk,
                        variables = c(region,Country,sex,female_lab,age,currsmok,bpms_hypt_med,sbp_avg_rev,bmicat_lab,edulte4yr,wealth_quintile,agegrp))

tableS2.d<-tableS2.svy %>% summarise(pctfem=survey_mean((female_lab=="Female" & !is.na(female_lab)),proportion=TRUE, prop_method="asin"),pctm=survey_mean((female_lab=="Male" & !is.na(female_lab)),proportion=TRUE, prop_method="asin"),pctsexmiss=survey_mean((is.na(female_lab)),proportion=TRUE, prop_method="asin"),pct30to34=survey_mean((agegrp=="30-34" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct35to44=survey_mean((agegrp=="35-44" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct45to54=survey_mean((agegrp=="45-54" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct55to64=survey_mean((agegrp=="55-64" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"), pct65plus=survey_mean((agegrp=="65+" & !is.na(agegrp)),proportion=TRUE, prop_method="asin"),pctagemiss=survey_mean(is.na(agegrp),proportion=TRUE, prop_method="asin"), pctunderwt=survey_mean((bmicat_lab=="<18.5" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"), pctnormal=survey_mean((bmicat_lab=="18.5-24.9" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"), pctoverwt=survey_mean((bmicat_lab=="25.0-29.9" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"),pctobese=survey_mean((bmicat_lab==">=30.0" & !is.na(bmicat_lab)),proportion=TRUE, prop_method="asin"),pctbmimiss=survey_mean(is.na(bmicat_lab),proportion=TRUE, prop_method="asin"),pctnonsmoker=survey_mean(currsmok==0 & !is.na(currsmok),proportion=TRUE, prop_method="asin"),pctsmoker=survey_mean(currsmok==1 & !is.na(currsmok),proportion=TRUE, prop_method="asin"),pctsmokemiss=survey_mean(is.na(currsmok),proportion=TRUE, prop_method="asin"),pctedultp=survey_mean(edulte4yr==1 & !is.na(edulte4yr),proportion=TRUE, prop_method="asin"),pcteduprimplus=survey_mean(edulte4yr==0 & !is.na(edulte4yr),proportion=TRUE, prop_method="asin"),pctedumiss=survey_mean(is.na(edulte4yr),proportion=TRUE, prop_method="asin"), pctWQ1=survey_mean(wealth_quintile==1 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQ2=survey_mean(wealth_quintile==2 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQ3=survey_mean(wealth_quintile==3 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"), pctWQ4=survey_mean(wealth_quintile==4 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQ5=survey_mean(wealth_quintile==5 & !is.na(wealth_quintile),proportion=TRUE, prop_method="asin"),pctWQmiss=survey_mean(is.na(wealth_quintile),proportion=TRUE, prop_method="asin")) %>% dplyr::select(pctfem,pctm,pctsexmiss,pct30to34,pct35to44,pct45to54,pct55to64,pct65plus,pctagemiss,pctunderwt,pctnormal,pctoverwt,pctobese,pctbmimiss,pctnonsmoker,pctsmoker,pctsmokemiss,pctedultp,pcteduprimplus,pctedumiss,pctWQ1,pctWQ2,pctWQ3,pctWQ4,pctWQ5,pctWQmiss) %>% mutate_all(funs( round(.* 100,1))) 
#t(bind_rows(tableS2.d,tableS2.m,tableS2.m.1))
#reported in manuscript
t(bind_rows(tableS2.m.1,tableS2.d))

dim(mdat);dim(dat)
```



```{r : Fig S5 - mult-female}
dat_female<-dat%>% filter(female_lab=="Female")
md<-dat_female %>% group_by(Country) %>% summarise_at(.vars= vars( edulte4yr,poor,wpartner,working), funs(miss=sum(is.na(.))))%>% left_join(.,dat_female %>% group_by(Country) %>% summarise(n=n(),n.complete=sum(!is.na(edulte4yr) & !is.na(poor) & !is.na(wpartner) & !is.na(working)  ))) 
#%>% mutate_at(.vars= vars( edulte4yr_miss, poor_miss,wpartner_miss,working_miss,n.complete),funs(paste(.," (",sprintf("%.1f",.*100/n),")",sep=""))) 
countrylist<-dat_female %>% group_by(region,Country) %>% summarise(n=n()) %>% pull(Country)
est<-data.frame(country=rep(countrylist,each=4),rf=rep(c("edulte4yr", "poor","wpartner","working"), length(countrylist)),setNames(data.frame(matrix(NA,nrow=4*length(countrylist),ncol=2)),c("beta","se")))
countrylist_rev<-countrylist[!(countrylist %in% c("Uganda","Indonesia","Kazakhstan",unique(with(md[md$n.complete==0,],paste(Country)))))]
i<-0
for (country in countrylist_rev){
  i<-grep(country,countrylist)
  print(i)
  anadat<-na.omit(dat_female[dat_female$Country == country,c("Country","WHOrisk","edulte4yr","poor","psu_num","weq_CVDrisk","wpartner","working")])
  if (dim(anadat)[1]>0){
    if (!(country %in% c("Uganda","Grenada","Morocco"))){
      if (sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0) {m<-lmer(log(WHOrisk)~I(edulte4yr==0)+I(poor==0)+wpartner+working+(1 | psu_num),data = anadat, weights = weq_CVDrisk)

if (isSingular(m, tol = 1e-05)) print(paste ("singular fit - F&M - ",country)) else est[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[-1,1:2]
    } 
        
  }
  else {
    if (sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0) {m<-lm(log(WHOrisk)~I(edulte4yr==0)+I(poor==0)+wpartner+working,data = anadat, weights = weq_CVDrisk)
est[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[-1,1:2]
    } 
  }
}}
estval<-data.frame(est,pctchange=round(100*(exp(est$beta)-1),1),pctchange_l=round(100*(exp(est$beta-1.96*est$se)-1),1),pctchange_u=round(100*(exp(est$beta+1.96*est$se)-1),1))

####   plot 
regcoef <- read.xlsx("TabsFigs.xlsx",sheet="multi_F",colNames=TRUE)
regcoef<-regcoef[!(regcoef$country %in% unique(regcoef$country[ regcoef$country%in% dat$Country & is.na(regcoef$beta_F)])),]
dim(regcoef)

## prim school or above
fpdat<-regcoef[regcoef$level=="edulte4yr" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,5),TRUE,rep(FALSE,7),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,12)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  middle/richer/richest
fpdat<-regcoef[regcoef$level=="poor" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,5),TRUE,rep(FALSE,7),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,12)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  married
fpdat<-regcoef[regcoef$level=="wpartner" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,5),TRUE,rep(FALSE,7),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,12)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

# working
fpdat<-regcoef[regcoef$level=="working" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
#xtcks<-c(-66,0,9)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,5),TRUE,rep(FALSE,7),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,12)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=c(-71,0,10),align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

```

```{r : Fig S6 - multi- male}
dat_male<-dat%>% filter(female_lab=="Male")
md<-dat_male %>% group_by(Country) %>% summarise_at(.vars= vars( edulte4yr,poor,wpartner,working), funs(miss=sum(is.na(.))))%>% left_join(.,dat_male %>% group_by(Country) %>% summarise(n=n(),n.complete=sum(!is.na(edulte4yr) & !is.na(poor) & !is.na(wpartner) & !is.na(working)  ))) 
#%>% mutate_at(.vars= vars( edulte4yr_miss, poor_miss,wpartner_miss,working_miss,n.complete),funs(paste(.," (",sprintf("%.1f",.*100/n),")",sep=""))) 
countrylist<-dat_male %>% group_by(region,Country) %>% summarise(n=n()) %>% pull(Country)
est<-data.frame(country=rep(countrylist,each=4),rf=rep(c("edulte4yr", "poor","wpartner","working"), length(countrylist)),setNames(data.frame(matrix(NA,nrow=4*length(countrylist),ncol=2)),c("beta","se")))
countrylist_rev<-countrylist[!(countrylist %in% c("Uganda","Indonesia","Tajikistan",unique(with(md[md$n.complete==0,],paste(Country)))))]
i<-0
for (country in countrylist_rev){
  i<-grep(country,countrylist)
  print(i)
  anadat<-na.omit(dat_male[dat_male$Country == country,c("Country","WHOrisk","edulte4yr","poor","psu_num","weq_CVDrisk","wpartner","working")])
  if (dim(anadat)[1]>0){
    if (!(country %in% c("Uganda","Grenada","Morocco"))){
      if (sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0) {m<-lmer(log(WHOrisk)~I(edulte4yr==0)+I(poor==0)+wpartner+working+(1 | psu_num),data = anadat, weights = weq_CVDrisk)

if (isSingular(m, tol = 1e-05)) print(paste ("singular fit - F&M - ",country)) else est[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[-1,1:2]
    } 
        
  }
  else {
    if (sum(table(anadat$edulte4yr)<5)==0 & sum(table(anadat$poor )<5)==0 & sum(table(anadat$wpartner )<5)==0 & sum(table(anadat$working )<5)==0) {m<-lm(log(WHOrisk)~I(edulte4yr==0)+I(poor==0)+wpartner+working,data = anadat, weights = weq_CVDrisk)
est[(4*i-3):(4*i),3:4]<-summary(m)$coefficients[-1,1:2]
    } 
  }
}}
estval<-data.frame(est,pctchange=round(100*(exp(est$beta)-1),1),pctchange_l=round(100*(exp(est$beta-1.96*est$se)-1),1),pctchange_u=round(100*(exp(est$beta+1.96*est$se)-1),1))

####   plot 
regcoef <- read.xlsx("TabsFigs.xlsx",sheet="multi_M",colNames=TRUE)
regcoef<-regcoef[!(regcoef$country %in% unique(regcoef$country[ regcoef$country%in% dat$Country & is.na(regcoef$beta_F)])),]
dim(regcoef)

## prim school or above
fpdat<-regcoef[regcoef$level=="edulte4yr" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,3),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,11)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  middle/richer/richest
fpdat<-regcoef[regcoef$level=="poor" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,3),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,11)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

##  married
fpdat<-regcoef[regcoef$level=="wpartner" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,3),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,11)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=xtcks,align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

# working
fpdat<-regcoef[regcoef$level=="working" |is.na(regcoef$level)|is.na(regcoef$country),c("country","pctchange_F","pctchange_l_F","pctchange_u_F")]
xtcks<-round(c(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],0,0.5*range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2]))
#xtcks<-c(-66,0,9)
fpdat$country[fpdat$country %in% unique(dat$Country)] <- paste("  ",fpdat$country[fpdat$country %in% unique(dat$Country)]) 
#fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(pctchange_F," (",pctchange_l_F,", ", pctchange_u_F,")")))
fpdat$effect<-with(fpdat,ifelse(is.na(pctchange_F),"",paste(sprintf("%.0f",pctchange_F),"%",sep="")))
par(mfrow=c(1,1),mai=c(0.1,0.1,0.1,0.1))
forestplot(labeltext=fpdat[,c("country","effect")], graph.pos=2, graphwidth = unit(2, "inches"),
           mean=fpdat$pctchange_F, 
           lower=fpdat$pctchange_l_F, upper=fpdat$pctchange_u_F,
           title="",
           xlab="Relative change (%)       ",
#hrzl_lines=list( "2" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),                 "10" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922"),"14" = gpar(lwd=40, lineend="butt", columns=1, col="#99999922")), 
           is.summary=c(rep(TRUE,2),rep(FALSE,6),TRUE,rep(FALSE,3),rep(TRUE,2),rep(FALSE,7),TRUE,rep(FALSE,11)), txt_gp=fpTxtGp(label=gpar(cex=0.9),ticks=gpar(cex=1),xlab=gpar(cex = 1),title=gpar(cex = 1),summary = list(
             gpar(fontface="bold",cex=0.9))), 
#            txt_gp=fpTxtGp(label=gpar(cex=1.2),ticks=gpar(cex=0.9),xlab=gpar(cex = 1.2)),
           col=fpColors(box="black", lines="gray70", zero = "gray50"),
           zero=0, cex=0.9, lineheight = "auto", boxsize=0.3, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.3,
#xticks=round(sort(c(seq(range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[1],range(c(fpdat$pctchange_l_F,fpdat$pctchange_u_F),na.rm=T)[2],length.out=5),0)))
xticks=c(-65,0,5),align = c("l", "r"),mar = unit(c(-1,1,-7,0), "mm"))

```

```{r : numbers in manuscript text}
table1.svy <- dat %>% 
            as_survey_design(stratum = stratum_num,
                        ids = psu_num, 
                        weights = weq_CVDrisk,
                        variables = c(region,Country,sex,female_lab,age,currsmok,bpms_hypt_med,sbp_avg_rev,sbp_adjtrt,bmi))

t1<-table1.svy %>% group_by(Country) %>% summarise(meanage=survey_mean(age),medage=survey_quantile(age,0.5),pctsmok=survey_mean(currsmok, proportion=TRUE, prop_method="asin"),bmi=survey_quantile(bmi,c(0.25,0.5,0.75)),pctF=survey_mean(sex, proportion=TRUE, prop_method="asin"),pctmed=survey_mean(bpms_hypt_med, proportion=TRUE, prop_method="asin"))

#age
round(c(quantile(t1$meanage),quantile(t1$meanage,0.75)-quantile(t1$meanage,0.25),range(t1$meanage,na.rm=T)),1)
#female
round(100*c(quantile(t1$pctF),quantile(t1$pctF,0.75)-quantile(t1$pctF,0.25),range(t1$pctF,na.rm=T)),1)
#female-unwtd
round(100*c(quantile(t1$pctF),quantile(t1$pctF,0.75)-quantile(t1$pctF,0.25),range(t1$pctF,na.rm=T)),1)

#smokers
round(100*c(quantile(t1$pctsmok),quantile(t1$pctsmok)[4]-quantile(t1$pctsmok)[2],range(t1$pctsmok,na.rm=T)),1)
#bmi
round(c(quantile(t1$bmi_q50),quantile(t1$bmi_q50)[4]-quantile(t1$bmi_q50)[2],range(t1$bmi_q50,na.rm=T)),1)
#med
range(t1$pctmed)
head(t1$Country[order(t1$pctmed)]);tail(t1$Country[order(t1$pctmed)])
round(100*c(quantile(t1$pctmed),quantile(t1$pctmed)[4]-quantile(t1$pctmed)[2],range(t1$pctmed,na.rm=T)),1)

t2<-table1.svy %>% subset(.,bpms_hypt_med==0) %>% group_by(Country) %>% summarise(sbp_notrt=survey_quantile(sbp_avg_rev,c(0.25,0.5,0.75)))
round(c(quantile(t2$sbp_notrt_q50),quantile(t2$sbp_notrt_q50)[4]-quantile(t2$sbp_notrt_q50)[2],range(t2$sbp_notrt_q50,na.rm=T)),1)

t3<-table1.svy %>% subset(.,bpms_hypt_med==1) %>% group_by(Country) %>% summarise(sbp_trt=survey_quantile(sbp_avg_rev,c(0.25,0.5,0.75)))
round(c(quantile(t3$sbp_trt_q50),quantile(t3$sbp_trt_q50)[4]-quantile(t3$sbp_trt_q50)[2],range(t3$sbp_trt_q50,na.rm=T)),1)

t4<-table1.svy %>% subset(.,bpms_hypt_med==1) %>% group_by(Country) %>% summarise(sbp_trt=survey_quantile(sbp_adjtrt,c(0.25,0.5,0.75)))
round(c(quantile(t4$sbp_trt_q50),quantile(t4$sbp_trt_q50)[4]-quantile(t4$sbp_trt_q50)[2],range(t4$sbp_trt_q50,na.rm=T)),1)


# med 10-yr CVD risk (IQR)
avCVDrisk_men
avCVDrisk_women


#percentage of people indicated for BP medication use who reported actual use of a BP medication 
avBPmed_indicatedMen
avBPmed_indicatedWomen

# proportion of people those taking BP medication were not indicated for such based on their CVD risk status
avNoIndicated

#median percentage of people overall indicated for BP medication who were taking a BP medication 
trt_svy <- dat %>% filter(is.na(bpms_hypt_med)==F & is.na(guideRx)==F ) %>% 
            as_survey_design(stratum = stratum_num,ids = psu_num, 
                        weights = weq_CVDrisk,variables = c(region,female_lab,Country,guideRx, 
                                    bpms_hypt_med))
table_trt0 <-trt_svy%>% 
  group_by(region,Country,guideRx) %>% 
  summarize(trtsvy = survey_mean(bpms_hypt_med,proportion=TRUE, prop_method="asin", 
                                vartype = "ci"))  %>% mutate(trtsvy=100*trtsvy,trtsvy_low=100*trtsvy_low,trtsvy_upp=100*trtsvy_upp)
quantile(table_trt0$trtsvy[table_trt0$guideRx=="indicated"])
round(c(27.997817,37.258871-16.290553 ),1)
```
